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FOREWORD 


This report presents four papers resulting from research 
conducted under a grant from NASA to the Old Dominion University 
Research Foundation entitled: "A Research Participation Program 

for Minority Engineering Students". The three undergraduate 
engineering students. Dale O. Douglas, Donna E. Holzmacher, 
and Zoa C. Lane, worked under the direction of Dr. Earl A.' 
Thornton, Associate Professor of Mechanical Engineering and 
Mechanics . 

Phe student— faculty team began their research in analysis 
of composite materials at Langley Research Center during a 
ten— week period in the summer of- 1974'. The work was continued 
during the academic year 1974—1975 at Old Dominion University. 

Dr. John G. Davis, Jr., - of the ' Composites Sectiony 
Materials Application Branch of the Materials Division served 
as technical monitor for- the program. For his cooperation, 
encouragement, and counsel the authors express their deepest 
appreciation. 
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FINITE ELEMENT ANALYSIS OF A PICTURE FRAME SHEAR TEST 

By 

Dale 0. Douglas 
INTRODUCTION 


Shear testing of composite materials is generally 'concerned 
with two principal areas of interest: (1) to determine the in- 

plane shear properties, or (2) to determine the interlaminar or 
normal shear properties. In-plane shear properties of a laminate 
are among the most difficult to determine because of problems in 
applying a. state of uniform shearing stress. Concepts for deter- 
mining in-plane shear properties include torsion tube tests, 
rail shear tests, and picture frame shear tests 

The most, direct method of applying pure shear is' by torsion 
of a tube. This test method has proven to be a reliable means 
of determining in-plane shear properties (ref. 1) . However 
fabrication techniques . for high quality ± 45° metal matrix 
composite tube’s have not yet been- established. ‘ The " difficulty 
of -fabricating high' quality tubes has stimulated research in 
other methods of shear testing. 

Another type of shear test is the rail shear test. It uses 
a thin laminate, loaded along its length by two pairs of rails, 
leaving -an unsupported central test section 

In the present study an analysis of a picture frame shear 
test performed at Langley Research Center is presented. The 
purposes of the study were to determine the stress distributions' 
in the picture frame shear test specimen and to determine the 
effect of local reinforcements on the stress distributions. 

DESCRIPTION OF TEST 


The experimental setup for a picture frame shear test is 
shown in figure 1. The picture frame shear test was used to 



produce in-plane shear stress in the test panel. The shear panel 
was bonded to a frame constructed from four 1 in.' x 1 in. steel 
edge bars designed to simulate fully clamped edge conditions. 

The panel specimen was bolted to a test frame by 0. 375-in. - 
diameter bolts , seven per side. At each corner of the test 
frame, loads were applied to the pin joints by the testing 
machine. Tensile loads were applied to the vertical pins, 
and compressive loads were applied to the horizontal pins to 
produce shear loading in the test specimen. 

TEST SPECIMEN 

The test specimens were made using 7 in. x 7 in. borsic 
aluminum sandwich shear panels. With the addition of 1 -in. x 
1 in. steel- edge bars, the overall dimensions of the shear panel 
specimen were 9 in. x 9 in. with a nominal thickness of 1- in. 

To permit installation of the pins on the test frame, a portion 
of the shear panel was cut away at each corner. Each corner had 
a radius of 0.25 in. The test specimen .is shown schematically in 
figure 2, 

The sandwich panel consisted of two face sheets separated by 
a- -honeycomb core. On each face sheet there were four plies 
(0.0285 in. thick) at a ± 45° layup. The panel face sheets were 
cut from 10-in. -square laminates. The filaments - of the laminate, 
were parallel to the applied loads. Some of the test specimens 
were reinforced with titanium doublers X 0.060 in. thick)' in the' • 
vicinity of the corner radii. 

ANALYSIS- OF SHEAR TEST 


Finite element analyses were made to determine the in-plane 
stress distributions in the shear panel. The finite models rep- 
resented the shear panel specimens using orthotropic, two-dimen- 
sional plane stress elements. Two general purpose finite element 
computer programs were utilized in the. analysis of the shear panel. 
The first was NASTRAN (NASA Structural Analysis Program) which was 
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used on the CDC-6600 computer at Langley Research Center. 

NASTRAN (ref. 2) is a general purpose digital computer programs, 
for the analysis of large complex structures. ' The second 
program, SAP IV (Structural Analysis Program) , was executed 
on an IBM-370, Model 158 computer at Virginia Polytechnic 
Institute & State University through the computer, center at 
Old Dominion University. SAP IV (ref. 3) is a structural 
analysis program for static and dynamic response of linear 
systems. Symmetry of loading, geometry, and material properties 
made the analysis of only- one quarter of the specimen sufficients 

NASTRAN embodied a finite element approach, wherein the 
distributed physical properties of the shear specimen were repre- 
sented by a model (fig. ■ 3) consisting of 490 membrane elements- 
that were interconnected at 529 grid points. * The grid point 
definition formed the basic framework for the structural models 
All other parts of the structural model were referenced either:- 
directly or indirectly to- the grid points.' Each grid point had,- 
two degrees of freedom, the in-plane displacements. The elementsr- 
used .in the analysis were the quadrilateral membrane element 
CQDMEM and the. triangular membrane element CTRMEM. 

- The steel edge bars of the test specimen were represented; . 
in "NASTRAN as rigid boundaries. The rigid boundaries were 
modeled using multipoint constraints in the NASTRAN program 
The constraint's were applied to grid points on the test frame 
edge of the finite element model so that these points deformed- 
as a straight line. Static loads' were applied to the structura] 
model through nodes constrained to the rigid boundary 

The loads were from Langley Research Center Test 560, -.Run 7 r 
a horizontal load of 5004.9 lb and a vertical load of 5039.4 lb ' 
are shown in figure 3 at the points of application. 

SAP embodied a finite element approach where the shear ' 
specimen was represented by a model (fig. 4) consisting of 554 
membrane elements that were interconnected at 595 nodal points- 
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The steel edge bars of the test specimen were represented in 
SAP as deformable boundaries. The deformable boundaries were 
simulated by the addition of 64 plane stress membrane elements 
to the NASTRAN model. The horizontal and vertical applied loads 
were represented by- .statically equivalent loads applied along the 
simulated boundary. Nine colinear loads were applied at nodal 
points nearest the center of each bolt hole. These loads were 
applied at an angle of 45 degrees. The magnitudes of these 
applied loads are given in figure 4. Stresses were computed at 
the centroid of each element using the stress print, option 
available in SAP. ’ 

The titanium doublers used for local reinforcement at corner 
radii were modeled with an addition of 21 finite elements' on 
existing elements at the extreme corner of the- sandwich panel. 

The material elasticity matrix for titanium and bor sic aluminum 
is given in table 1. 


fable 1. • Material elasticity matrix. 
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The NASTRAN finite element model of the shear panel, simu- 
lating a rigid boundary, had 1000 degrees of freedom. Using a 
CDC-6400 computer, it took 945 CPU seconds for the program to 
execute. In contrast to the NASTRAN model, the SAP finite model ' 



had 1132 degrees of freedom with a bandwidth of 1106. Due to 
the excessive storage required by the large bandwidth, the SAP- 
finite element program. was unable to execute. To optimize the- 
bandwidth, the nodes were renumbered using a computer program ,. 
BANSAP. With this renumbering, the SAP program had a final 
bandwidth of 69. It. then completed execution in 160 CPU seconds:^. 

RESULTS AND DISCUSSION 

The stresses computed in the shear panel for the loads 
^■Pplisd to the rigid (NASTRAN) and deformable (SAP) boundary 
models are given in figures 5 through 8. Normal stresses' a x 
and o are plotted as ordinates with the horizontal and 
vertical coordinates from the center of the shear panel as 
abscissae. 

The stress distributions along the horizontal ' and vertical. . 
axes of both the rigid and deformable boundary models are uniform? 
near the center of the specimen. The uniform stress values differs- 
considerably between. the two models; the uniform normal stpesr~~~ 
predicted -‘by the rigid boundary model, are nearly three times 
the stresses predicted" as suming a deformable boundary. These- 
results indicate that the assumption of a rigid boundary shouIcL- 
not be made. 

There is an appreciable stress concentration at the cornerr 
fillets. . The stress components perpendicular to .the lines of 
symmetry rise sharply at ‘the corners. For example, figure 7 
shows that : in the deformable boundary model the stress component 
cr y • increases from 10 ksi to about 105 ksi indicating a stress 
concentration factor of over 10. 

Contour plots of the principal shearing stresses for the 
rigid and deformable boundary models are shown in figures 9 and. 

10. The shearing stresses are uniform only oyer a small portion 
of the specimen. Figure 10 shows that the shearing stress may 
vary by as much as 25 percent over the center one-half of the 
specimen. 



The effect' of the reinforcing titanium doubler on the normal 
stresses and <x is shown in figures 11 through 14. These 

results indicate that the doubler significantly reduced the 
stresses along the x-axis near the fillet. The critical stress 
Oy oft this axis was reduced by about one-half. However, stress 
distributions along the vertical axis and in the center of the 
shear panel show no reinforcing effects of the titanium doubler. 

The contours of the principal shearing stress in the specimen 
with the titanium doubler are shown in figure 15. By comparing 
this figure with figure 10 it can be seen that the doubler tended 
■to reduce the region of nearly uniform shearing stress since the 
contours in figure 15 are closer to the center of the panel. As 
expected there is also an appreciable local disturbance in the 
shearing stresses in the- vicinity of the doubler. 

CONCLUDING REMARKS 


Two finite element analyses of a picture ' frame shear test 
of a borsic. aluminum test specimen have been performed. Two. 
methods for modeling the. specimen test frame have been investi- 
gated. ~ Results for nominal stresses and principal shear stress 
have been presented for Test- 560, Run 7 conducted at Langley 
Research Center. 

~ There were striking differences in the stress distributions 
predicted by -the rigid (NASTRAN) and ' deformable (SAP) boundary 
models.. It was found that it is’ not realistic to assume the : 
test fixture to be a. rigid frame. In the , regions - of nearly 
uniform stress, the stresses predicted by the deformable- 
boundary models" were approximately one third of the stresses 
predicted by the rigid boundary model. In the vicinity of the 
corner, the stresses predicted by the two models nearly coincided. 

The constant principal shear stress, t was uniform 

incix 

over only a very small region in the center of the shear panel 
specimen. Moreover, at the corners near the fillets, there were 
steep gradients with stresses being highly concentrated. 
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The effect of a local reinforcing titanium doubler has 
been evaluated. It was found that the doubler reduced the 
maximum nominal stress in. the vicinity of the fillet by about 
50 percent. 
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Figure 1. Picture Frame Shear Test Experimental Setup at Langley Research Center, 
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Figure 2. Schematic of Test Specimen. 
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Figure 5. Normal Stress a as a Function of Horizontal 
Coordinate Along x Center Line of Shear Panel 
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Ficrure 10. Contours of Constant Principal Shear 
Stress, T max ' Predicted by SAP . 

(Deformable Boundary) . 
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Figure 12 . Normal Stress a as a Function of Vertical 
Coordinate Along x Center Line of Shear Panel 
Specimen. 
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the centroidal coordinate of the finite elements bordering the 
hole. Comparison of figures 6 through 8 and figure 10 shows 
that qualitatively the finite element analysis of the anistropic 
composite and the isotropic elasticity solution are in close 
agreement. This agreement serves to validate the finite element, 
solution. 

The variation or the longitudinal membrane force in an iso- 
tropic infinite medium is shown in figure 11 in terms of the x 
coordinate of the composite specimen to facilitate comparison 
with the finite element solution given in figure 9. The elas- 
ticity solution shows an extremely sharp gradient for the mem- 
brane force in the vicinity of the hole. This sharp variation 
raises questions about the accuracy of the finite element solution 
in this region. Since the NASTRAN finite element assumed constant 
stress within the element, it is possible that the peak stress" 
was underestimated because not enough elements were used 'to", 
accurately represent the stress gradient. The variation of -the 
stress away from the hole according to .the isotropic solution 
shows that in a distance of about five" radii (5a = 0 ; 48;~in. ) 
away from the -hole the force has decreased to one-tenth o.f its 
maximum value. This result supports the findings -of figures 
6 through 8 in which the membrane force distributions in the’ 
center and outside holes were very nearly the same. This 
occurred because there were no hole interaction effects since 
the holes were more than five radii apart. Only very small 
edge effects .were present for the same reason. 

CONCLUDING REMARKS 

A finite element analysis of an extra graph 

bolted joint specimen has been performed. Two methods were used 
to represent bolt transfer loads. The first method assumed a 
perfect fit and modeled the bolt loading as a cosine distributchn 
over one-half of the boundary of the hole. The second method 
assumed an imperfect fit and used a nonlinear computer analysis 
to determine the contact area and bolt transfer loads. The 
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Figure 13. Normal Stress a as a Function of Horizontal 
Coordinate Along y Center Line of Shear Panel 
Specimen . 
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Figure 15 . Contours of Constant Principal Shear Stress , 
T max / Predicted by SAP (Deformable Boundary 

with Reinforcing Titanium Doubler) . 
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BANSAP: A BANDWIDTH REDUCTION PROGRAM FOR SAP IV 

' By 

Donna E. Holzmacher 
INTRODUCTION 

For analysis, a structure may be broken down into parts 
known as finite elements. The elements of the structure may be 
one-dimensional such as a rod, two-dimensional such as a triangle 
or quadrilateral, or three-dimensional such as a parallelepiped. 
The elements are positioned and described by nodes which, when - 
connected, describe the structure. Static analysis using finite' 
elements^ is- accomplished by solving simultaneous equations. 

These equations when written in matrix form are characterized 
by banded coefficient matrices. Computer time and storage- car 
be saved if the bandwidth of the matrix is a minimum. This 
occurs with adept "numbering of the nodal points of' the ..structure. 
If the nodes are numbered in an optimum way the non-zero values . - 
in the matrix will lie in a band about the diagonal. The band- 
width of a' matrix is defined here as -the maximum difference 
between any two connected nodes plus one to take into account 
the .diagonal term. 

As a particular example, consider the plane" structure shown 
in" figure 1. The displacements, of this structure are "determined 
by solving 

(K) {U} = {P} 

where (K) is the stiffness matrix, {U> is the displacement 
vector, and {P} is the load vector. The (K) matrix is 
arranged according to the connectivity of the nodes 1 through 9 
of the triangular elements. The connectivity matrix for the 
above structure is represented in figure 2 showing that node 1 
is connected to nodes 2, 8, and 9, and node 2 is connected to 



nodes 1, 2, 3, 4, and 9, etc. The actual values in the stiffness 
matrix, corresponding to the positions of the matrix above, 
depend on the geometry and material of the structure. 

The bandwidth of the connectivity matrix shown in figure 2 
is 9. If the nodes are renumbered as in figure 3, the corres- 
ponding connectivity matrix as shown in figure 4 has a new, 
reduced bandwidth of 4. 

In order to efficiently renumber the nodes of structures 
for finite element analysis a number of algorithms have been 
developed and incorporated into bandwidth reduction programs . 

Prior to 1969, authors who developed techniques to reduce the 
bandwidths of matrices included Always and Martin, Tewarson, 

Rosen, and Akyus and Utku (refs. 1 through 4>. In 1969, 

Cuthill and McKee's (ref. 5) algorithm arranged the rows of -the 
connectivity matrix with regard to the increasing number 'of non- 
zero off-diagonal elements. This algorithm was used in a program 
called BANDIT which serves a a preprocessor for NASTRAN. 

H.R. Groom's algorithm for bandwidth reduction was intro- 
duced in 1972 (ref. 6).. Groom systematically moved closer to- 
gether rows and columns which were far apart and coupled." 

In 1973, Collins (ref. 7) presented the algorithm upon which 
the program,. ’BANSAP , developed in this study is based. After 
work on BANSAP had begun, Rodrigues' (ref. 8) presented a new’ 
algorithm which, for two sample problems presented, showed a 
smaller bandwidth than the Cuthill and McKee, the Groom, or .the 
Collins algorithms. 

The objective of this paper is to describe a study- undertaken 
to incorporate the Collins bandwidth algorithm in a data pre- 
processing computer program for the finite element program SAP IV 
(Structural Analysis Program - IV) . First to be presented will 
be Collins ' algorithm for bandwidth reduction which contains two 
subroutines, SETUP and OPTNUM. A description of the SAP IV 
preprocessing program BANSAP will then be given including its 
capabilities and limitations. Results from application of the 
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program to example problems will be presented and discussed 
User instructions, the BANSAP program listing, and sample 
output are presented in. appendices. 


COLLINS BANDWIDTH REDUCTION ALGORITHM 

Collins’ algorithm for bandwidth reduction includes two 
subroutines, SETUP and OPTNUM. His procedure shall be illus- 
trated using the structure in figure 1. 

In the first subroutine, SETUP, a list is generated showing 
the connections between the different nodes shown in figure 1. 

The relations established by SETUP are displayed in table 1. 

The information is stored in arrays suitable for use in sub- 
routine OPTNUM. - The subroutine SETUP also determines the 
original bandwidth of the structure. 

The subroutine used to renumber related nodes is OPTNUM. 
OPTNUM locates the origin of the different numbering schemes 
at each node in turn, making' the number of permutations of ' 
schemes equal- to the number of' nodes. In other words, OPTNUM 
first renumbers the nodes around old node number one making old 
node number one the origin of the new scheme. / OPTNUM then deter- 
mines the bandwidth of this scheme. Next OPTNUM goes to old 
node number two and starts its new origin in the position of 
this node. It renumbers the nodes connecting node' two, one at 
a time, and determines the. maximum difference between the new 
connected nodes. If the maximum difference is less than the - 
lowest maximum difference of the preceding schemes, it continues- 
with the renumbering until the scheme is complete. If not, the 
current scheme is abandoned. After completion or abandonment • 
of a scheme, OPTNUM proceeds to the next scheme starting with 
a new origin at the next sequential old node number. The 
scheme which is retained by OPTNUM is that which exhibits the 
lowest maximum difference between related nodes. The sequence 
of renumbering schemes for figure 1 is shown in table 2. 



Collins' algorithm is set up to handle the renumbering of 
nodes for elements containing up to four nodes. Reference 5 
indicates that this method has been applied to solid elements 
but not very successfully. 


SAP IV PREPROCESSING PROGRAM 

A program, BANSAP, has been written using the Collins 
algorithm as a preprocessing program for SAP IV. BANSAP con- 
sists of four subroutines: SAPIN, SETUP, OPTNUM, and SAPOUT 

as shown in figure 5. 

The first subroutine, SAPIN, reads the data in the formats 
stipulated by SAP IV and stores element and node connections 
according. to type. BANSAP. is set up to handle two basic types 
.of finite elements : elements connecting two nodes, and elements 

connecting three or four nodes. The two node elements which can 
be entered into subroutine SAPIN- are either the truss, beam, or' 
boundary. The actual renumbering of a two node element is the 
same for either element. The only difference in the. handling 
of these elements by BANSAP is in their SAP IV formats.”. The 
three or four node elements which may be entered into subroutine 
SAPIN include membranes, axisymmetric two-dimensional elements 
and plate bending elements. ' Again,' the only difference in the 
handling of these three and four node elements is in their 
SAP IV formats . If more than one type' of element comprises the 
structure, the elements may be grouped according to their type. 
As' is required for SAP IV, nodes must, be sequentially numbered 
from one . ' 

From subroutine SAPIN, BANSAP goes on to subroutines SETUP 
and OPTNUM. The new bandwidth is printed and a list of old 
number node numbers and new numbers is generated. As a user 
option the subroutine SAPOUT will punch the original elements 
with the new node numbers. Program BANSAP has been dimensioned 
in this paper to permit up to 1000 nodes and 1000 elements. 



APPLICATIONS OF BANSAP 


Applications of BANSAP are presented in table 3. The first 
two problems shown are illustrations of reduction in bandwidth 
which may be attained for simple problems. Problem 3 is an 
example taken from structural analysis of a ship radar tower. 

The last two entries are practical problems encountered in 
finite analysis of composite material structures. 

The first illustration is the sample finite element scheme 
shown in figure 1. After renumbering by BANSAP the bandwidth 
was reduced from 9 to 4 and the final scheme is shown in figure 3. 

The truss problem shown in figure 6a is a wagonwheel. Afte->~ 
processing by BANSAP the bandwidth was reduced from 9 to 6. ’It 
has been found, however, that this value is not the optimum' 
-bandwidth. Collins has noted that the wagonwheel problem is 
a. special case and the true optimum bandwidth occurs' when the 
node -number of the hub of the wheel is set equal to half the 
number of spokes plus one. The optimum bandwidth of the wagon- 
wheel shown in figure '6 is actually 5. 

The third structure is the ship's radar tower shown in 
figure 7 . ' The original numbering scheme shown is nearly optimum’ 
with a bandwidth of 12 since the renumbering scheme only reduces • 
the bandwidth to 9. For such structures there is no appreciable 
gain by using BANSAP as -the ‘structures could- easily be numbered 
by -hand to obtain a small bandwidth.' 

The shear panel of figure '8 is an example of a greatly 
enlarged bandwidth which can occur from the addition of. new 
finite elements after the original .structure has been numbered. 
With the addition of new elements for the shear panel a bandwidth 
of 406 was obtained, but after BANSAP, the bandwidth was reduced 
to 35. 

The bolted joint specimen illustrated in figure 9 is a- good 
example of how BANSAP can be used to obtain an optimum bandwidth 
when the numbering scheme is difficult to select, by hand. The 
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nodes of the bolted joint specimen were originally numbered to 
permit easier data generation using a FORTRAN program. After 
the cards had 'been generated, BANSAP rerminbered the nodes to 
reduce the bandwidth from 168 to 28 . 

CONCLUDING REMARKS 

A FORTRAN program has been written for bandwidth reduction 
by nodal renumbering. The program is based upon the Collins 
algorithm and serves as a data preprocessor for the finite 
element program SAP IV. Applications of the preprocessing pro- 
gram to a number of simple and realistic problems have been 
presented. 

■Nodal renumbering for finite element .analysis may be required 
for a variety of reasons . Renumbering may be needed if new elements 
were to be added onto a previously numbered structure or if a. struc- 
ture is difficult to optimally number by hand. It may also be 
needed if the element and. nodal data were prepared by data genera- 
tion programs. Such reasons clearly show a need and use for a 
program such as BANSAP. 

BANSAP is an effective preprocessing program for SAP IV. 

The - algorithm used greatly reduces the bandwidth for reduced 
computer time and storage during the finite element analysis. - 
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APPENDIX A 
USER INSTRUCTIONS 
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USER INSTRUCTIONS 


CONTROL CARD (315) 

Columns 1-5 Number of different groups of elements 
6-10 Total number of nodes 
11 ” 15 The number zero for nopunched output and 
any number greater than zero for punched 
output 

The following types of elements are permitted in the program. 

Type 1 TRUSS 

CONTROL CARD (315) 

Columns 1-5 The number 1 

6-10 The number of e±emcm.a m yiuup ± 

Element Data Cards (315, 2A10) 

Columns 1-5 Element number 

6-10 Node number I 
11 - 15 .Node number J 

Type 2 BEAM 

CONTROL CARD (315) 

Columns 1-5 The number 2 

6-10 The number of elements, in group 2 

Element Data- Cards (-415,' 5A10) 

Columns .'1-5 Element number 
6-10 Node number - I 
11 - 15 . Node number J . 

16 - 20 Node number K; K' is any nodal point which 
■lies in the local i - 2> plane but not on 
the 1 axis (see ref. 9,. page iv.2.2) 

Type 3 MEMBRANE 

■ CONTROL CARD (315) 

Columns 1-5 The number 3 

6-10 The number of elements in group 3 

Element Data Cards (515, 5A10) 

Columns 1-5 Element number 

6-10 Node number I 

11 - 15 Node number J 

16 - 20 Node number K 

21 - 25 Node number L 
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Type 4 TWO D 


CONTROL CARD (315) 

Columns 1-5 The number 4 

6-10 The number of elements in group 4 

Element Data Cards (515, 5A10) 


Columns 


1 - 

5 

Element number 

6 - 

10 

Node 

number 

I 

11:- 

15 

Node 

number 

J 

16 - 

20 

Node 

number 

K 

21 •- 

25 

Node 

number 

L 


Type 6 PLATE 

CONTROL CARD (315) 

Columns 1-5 . The number 6 

6-10 Number of plate elements 

Element Data Cards (515, 5A10) 


Columns 


.1-5 


Element number 


6 - 

10 

Node 

-number 

I 

11 - 

15 

"Node 

number 

J 

16 - 

20 

Node 

number 

K 

21 - 

25 

Node 

number 

L 


Type 7 - BOUNDARY (LINEAR SPRING) 
CONTROL CARD (315) 


Columns 1-5 The number 7 ■ 

.6 - 10 The number of elements in group -7 . 

Element Data Cards (215, 6A10) 


Columns 1 
6 


5 Node N, at which the element is placed 
10 Node I 
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APPENDIX B 


BANSAP SOURCE LISTING 



000003 

000003 

000003 

000003 

000003 

000003 

000003 

000004 

000005 


000006 

000012 


000012 
3000 16 
0000 16 
0000 20 
000023 
0000 24 
000027 
000030 


PROGRAM 8ANSAP! I NPUT , OUTPUT , PUNCH, 

* TAP E5= INPUT, TAP E6=0UTPUT,TAPE7=PUNCH) 

DIMENSION NEW JT { 1000 )» JOINT ( 1000) - 

COMMON LMENTS, JT ( 4000) , MEMJT ( 3000) ,JMEM(1000) ,JNT(1000) 

COMMON/BAND/ 10 IFF, M INMAX 

COM.MON/CONTR/NELG , ITYPEl 5) , NE L( 5) , NODES, IPUNCH 
COMMON/JUNK/AI 1000,6) 

COMMON/UN IT/ IN, IT, IP 
IN = 5 
IT a 6 
IP = 7 
C 

C ' JMFM(T) = NUMBER OF NODES TO WHICH A SINGLE NODE IS CONNECTEC 

C JT(.I) = WORKING ARRAY 

C 'MEMJT(I) = IDENTITIES OF NODES TO WHICH A NODE IS CONNECTED 

C 

WRITECIT, 12) 


12 ! 
1 

FORMAT! IH 1 , 9 C /) , 
36X, 52HBBBBB 

AAAAA 

N 

N 

SSSSS 

AAAAA 


PPPPP/ 

2 

36X., 53HB B. 

A 

A 

NN 

N 

s 

A . 

A 

• P 

Pi 

3 

36X.53HB B 

A 

A • 

N 

N N 

s 

- A 

A 

P: 

Pi 

4 

36X, 52HBBBBB 

-AAAAA A A 

N 

N N' 

'*• SSSS 

AAAAAA A 

. PPPPP/ . 

5 

36X.52H3 B 

A 

A ‘ 

N 

N N 

s 

A 

A 

P 

/ 

- 6 

36X , 52H3 B.- 

A 

A 

N 

NN 

s 

A ' 

A 

P 

/ 

7 

36X, 52H3B88B 

A 

. A 

N 

N 

sssss 

A 

A 

P , 


WRITP (IT, 16) - • ' ’ 

16 F0RMAT(1H1,5X,19H INPUT DA T "A ,//,) 

" DO 10 1=1 ,1000 
10 JNT I I ) = 0 

DO 20 1=1,4000 
20 JT(I) = 0 

DO 30 1=1,3000 
30 MEMJT ( I ) =0 
C 

CALL SAPI-N 
C 

C SUBROUTINES SETUP AND OPTNUM FROM/ 

C -BANDWIDTH REDUCTION BY AUTOMATIC RENUMBERING-, R.J. COLLINS', 

C INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING 

C VOLUME 6, 1973, PP 345-356. 

C 


000034 


CALL SETJP 



000035 

C 


WRIT? [IT, 321 

000041 



DO 40 1= 1, NODES 

000043 



NO = JMFMU) 

000045 



L 1 = 8* ! I— 1 ) + 1 

000047 



L2 = LI *N0 -1 

000051 



WRITE! IT, 34) I , NO , l MEMJT! L) ,L=L1,L2> 

000070 


40 

CONTINUE 

000073 



MINMAX = IDTFF + 1 

0000 75 

r 


WRITE! IT, 36) MINMAX 

000102 

V 

r 


CALL OPTNUM 

000103 

U 


MINMAX = MINMAX +1 

000105 



WRITE (IT, 38) MINMAX 

000112 



WRITE! IT, 42) 

000116 

L» 

“ c 


CALL SAP3UT 

0001 17 

U 

32 

FORMAT I 1H1 , 12 X , 4HN0DE , AX, JM tM , 16X, SHMfcM J 1 ,/VJ 

000117 


34 

F0RMAT(11X,2I5,10X,9I6) 

0001 17 


36 

FORMAT!// , 20X, 20H0R I GI NAL BANDWIDTH = >14 1 

000117 


3 8 

FORMAT!// ,20X,14HNEW 3 ANDW IDTH= » 14) “ 

0001 17 


42 

FORMAT! 1HI, I0X,33H0LD NOOE NUMBER NEW NODE NUMBER,/) 

000117 



STOP 

000121 



END 



000002 

000002 

000002 

000002 

000002 


000002 

000014 

000026 

000030' 

000037 

000051 

000055 

000057 

000061 


C 

C 

C- 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


SUBROUTINE SAPIN 

COMMON LMENTS, JT { 4000 ), ME MJT { 8000 ) , JMEMli.000) , JNT( 1000) 

C 0MM0N/8AND/ I D I FF r M INMAX 

COM MON/CO NT R/ NELG , I TYPE ( 5 ) ,NELl 5 ) , NODES, IPUNCH 

COMMON/ JUNK /A 1 1000,6) 

COMMON/UNIT/ IN, IT, IP 


NELG = NUMBER OF DIFFERENT GROUPS OF ELEMENTS (LES S THAN 5 
NODES = TOTAL NUMBER OF NODES 

IPUNCH = ZERO FOR NO PUNCHED OUTPUT, NUMBER GREATER. THAN 
ZERO FOR PUNCHED OUTPUT 
N = ELEMENT TYPE 

NE‘ = NUMBER OF ELEMENTS OF TYPE N 

LMENT = TOTAL NUMBER OF ELEMENTS 

I TYPE { I ) = TYPE OF ELEMENT 

NEL(I) = NUMBER OF ELEMENTS IN A GROUP 

ITYPE ' ELEMENT NUMBER OF NODES 


1 TRUSS 2 

2 BEAM 2 

3 MEMBRANE 3 OR 4 

4 TWO 0 3 OR 4 

5 ’ BRICK 8 

6 PLATE 4 

7 BOUNORY 2 


READ .ELEMENT CARDS AND. STORE CONNECTIONS, 

R E A D ( I N ,"l 2 ) NELG, NODES, IPUNCH 

WR ITE ( I T, 14 ) .NELG, NODES, IPUNCH 

DO 200 11= 1 , NELG 

R EAD( IN, 1 2) N , NE 

WR I TE ( I T» 10 ) 1 1 , NE, N 

WR ITE ( IT, 50 ) 

rTYPEIII) =N 

NELUI) =NE 

LMENT S =0 

READ ELEMENT CONNECTIONS. FOR TRUSS, BEAM, OR BOUNOARV 
ELEMENTS ONLY TWO CONNECTIONS I AND J ARE NEEDED. FOR 



000062 

000063 

000076 

000076 

000116 

000116 

000117 
3001 17 
000137 
000137 

000140 

000140 

000164 

000164 

000165 

000165 

000171 

000171 

000172 

000172 

000212 

000212 

000212 

000214 

000216 

000220 

000221 

000223 

000226 

000230 

000232 


C ALL OTHER TYPES FOUR CONNECT EONS ARE POSSIBLE- I,J,K,L. ' 

C 

C STORE N03-E CONNECTIONS ACCORDING TO TYPEa 

C 

DO 210 JJ = 1 » NE 
GO TO ( 1, 2,3, 3, 5,3,7), N 
C 

1 CONTINUE 

READ ( I N >1 02') I,J,I AtJJ.U ,L=1,2) 

102 F0RMAT(5X,2I5,2A10) 

GO TO 300 
C 

2 CONTINUE 

READ! ! N ? 1 04 ) I,J,{ A { J J ,L) , L=1 » 6) 

104 F0RMATI5X ,2 13 , A5 , 5A10J 
GO TO 300 
C 

3 CONTINUE 

READ ( IN, 1 06 ) I,J,K,L,( AI JJ ,L l , L= l, 5) 

106 FORMAT { 5X ,415 * 5A10) 

GO TO 300 ' 

C 

5 CONTINUE 

WR ITEIIT, 103) 

108' ^OSMATIS'X .42HTHREE DIMENSIONAL ELEMENTS NOT IMPLEMENTED*/) 

GO TO 200 

e 

7 CONTINUE 

READ! IN»l 10 ) I,J, (A(JJ,L)»L=1,6) 

_ 110 F0RMA.T(2I5,6A10) 

300 CONTINUE 

IFCN.E.Q.7) GO TO 20 
III = JJ + LMENTS 

JJJ = III *■ 1000 ‘ 

JTCIII) = I 
JTCJJJ) = J 

I Ft N.LE. 2) GO TO 205 
K.KK = III + 2000 

LLL = III- + 3000 

JT(KKK) = K 
C 

C FOR TRIANGULAR ELEMENT SET REPEATED NODE NUM3ER EQUAL TO ZERO- 



000234 IFtK.EQ.U U=0 

300235 JT(LLL) = L 

000240 205 CONTINUE 

000240 WR I TE ( I T , 30 ) J J , J Tt I TI) , JT t JJ J) , JT ( XKK) , JT! LLL) 

000256 210 CONTINUE 

000261 LMENTS = LMENTS + NE 

000262 200 CONTINUE 

000265 50 FORMAT! 10 X , 7HELEMENT ,5X , 1H I ,7X, 1HJ,7X, 1HK,7X, 1HL,// J 

000265 10 FORMAT!/, 5X ,13HELEMENT GROUP, 12, 4H HAS,I3,17H ELEMENTS OF TYPE » 

l/) 

000265 12 FORMAT (315) 

000265 14 FORMAT! , 2 1 / J 

1 10X i 25HNUMBER OF ELEMENT TYPES =,I5,/ 

2 10X’, 25HNUM 8ER OF NODAL POLNTS = ,I5, / 

3- 10X, 25HPUNCHE0 ELEMENT CAROS =,I5»/ 

4 10X, 25H .EQ. 0 NO ,/ 

5 L0X.25H .EQ. 1 YES ,) - 

000265 30 FORMAT! 1QX, 15,418) 

000265 RETURN 

000266 END 
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030002 

000002 

000002 


000002 
000003 
000005 
0000 it 
000012 


000013 


0000 17 
000021 
000023 
000025 


C 

c 

c 

c 

c 

c 

c 

c 

c 

-C 


c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE SETUP 

COMMON LMENTS, JT ( 4000 }, MEMJT ( 8000 ) , JMEM(IOOO) ,JNTC10001 

COMMON/BA NO/I D IFF , MI NMAX 

COMMON /CQNTR/NELG, I TYPE ( 5 ) ,NELI 5) , NODES , I PUNCH 

NODES = TOTAL NUMBER OF NODES 

J NT I = ELEMENT NODE UNDER CONSIDERATION _ 

J SUB = LOCATION IN MEMJT(I) OF BEGINNING OF LIST OF 

NODES RELATED TO JNTI 
LMENT = TOTAL NUMBER OF ELEMENTS 

JMEM(I) -= NUMBER OF NODES TO WHICH A SINGLE NODE IS CONNECTS! 
MEMJTU) = IDENTITIES OF NODES TO WHICH A NODE I-S CONNECTED 
IDIFF. ■ = BANDWIDTH = IDIFF+1 FOR ORIGINAL SCHEME 

■ IDIFF = O' 

. DO 10 J= 1, NODES 
10 J MEM ( J ) = 0 

DO 60 J= 1, LMENTS 
DO 50 1= 1»’ 4 

NEXT STATEMENT DEPENDS ON THE NUMBER OF NODES FOR WHICH THE 
-PROGRAM IS DIMENSIONED. CURRENTLY' THE MAXIMUM NUMBER OF NODES 
IS 1000. 

JNTI = JTllOOO* (1-1) + J) 

IF JNTI EQUALS ZERO ALL NODES OF ELEMENT J HAVE 'BEEN 
CONSIDERED. 

IF( JNT.EQ.O) GO TO 60 
JSU3 = (JNTI - l) * 8' 

DO 40 1 1 = 1,4 •- ' . 

IF! II .EQ. I) GO TO 40 


C NEXT STATEMENT DEPENDS ON THE NUMBER OF NODES FDR WHICH THE 

C PROGRAM IS DIMENSIONED. CURRENTLY THE MAXIMUM NUMBER OF NOOES - 

C IS 1000. - 

C' 

C RELATED NODES ARE IDENTIFIED BELOW. 

C 


000027 

000033 


JJT = JT (1000 * t 1 1— X > +J ) 
IFCJJT.EQ.O) GO TO 50 



000035 

000036 

000040 

000041 


000044 
000046 
000050 
000052 
000055 
000062 
000064 
000066 
000071 
0000 72 


C DETERMINE WHETHER RELATIONSHIP BETWEEN JNTI AND JJT 

C HAS BEEN ESTABLISHED. 

c 

MEM1 = JMEMIJNTI) 

I F { MEM1 .EQ.O) GO TO 30 
DO 20 III = 1, MEM1 

I F C MEMJTl JSUB +III).EQ.JJT) GO TO 40 
C 

C FINO WIDTH OF ORIGINAL MATRIX BAND, 

C 

20 CONTINUE 

30 JMEMIJNTI) =JMEM( JNTI ) + 1 
IDUM JSUB + J M c M ( J NT E ) 

MEMJT{ IDUM) = JJT 

IF( IABSIJNT I-JJT) .GT.IDIFF) IDIFF = IABSIJNTI -JJT) 

40 CONTINUE 
50 CONTINUE 
60 CONTINUE - 
RETURN 
END 



000002 

000002 

000002 

000002 

000002 


000002 


000004- 

000005 


000006 

000007 


000013 
000013 
000014- 
000016 
0000 17 
000021 
000021 
000023 
000025 


000026 


SUBROUTINE OPT MUM 

DIMENSION NEWJT ( 1000 ) ? JOINTtlOOO) 

COMMON LMENTSt JT (4000 ), MEMJT (8000) , JMEMUOOQ) ,JNT!1000) 

COMMON/ BAND/ ID I FF,M INMAX 

COMMON/CONT R/NELG , ITYPE15) ,NEL( 5) , NODES, IPUNCH 
COMMON/UNIT/ INtIT, IP 

C JOINTII) = WORKING ARRAY 

C NEWJT (I) = WORKING ARRAY 

C ' JNT(I) = NEW NUMBERING SCHEME 

C — M INMAX = BANDWIDTH = MINHAX+1 FOR NEW SCHEME 

C 

C MINMAX IS INITIALIZED. - 

C - 

MINMAX = ID IFF. 

C 

C NEW SCHEME STARTS AT NODE OF OLD NODE NUMBER IK.- 

C 

DO 60 IK=1, NODES 
DO 20 J =1 , NODES 
C 

C JOINT(J) -AND NEWJT I.JJ INITIALIZED TO ZERO FOR EACH NEW 

C NUMBERING SCHEME. 

C 

JOINT! j)=- D 
20 NEWJT ! J ) = 0 
C 

C INITIAL IZE' FOR NEW NODE NUMBER ONE. 

C 

MAX =0 
1 = 1 

NEWJT(l) =IK 
JOt NT ( I KJ •=■ 1 
K = 1 

30 CONTINUE 

JDUM = NEWJT! II 
K4- = JMEMIJDUM) 

IF{K4.E0.0) GO TO 45 
C 

C LOCATE RELATED NODES IN MEMJTtlJ. 

C 

JSUB = (NEWJT! I ) -1) *8 


43 



000031 
000032 
000034 
000040 
000041 
0000 42 


000043 


000045 
000050 
000053 
000056 
000060 
000062 
000062 
000064 
000065- 
0000 71' 
000074 
000075 


DO- 40 JJ= 1,K4 

K5 = MEMJ T{ JSUB +JJ) 

IF{ J 0 1 NT ( K 5 ) .GT. 0 ) GO TO 40 
K = K + l 
NEW JT C K I =K5 
J0INTCK5}=K 
C 

C CHECK DIFFERENCE BETWEEN NEW NUMBERS OF RELATEO NODES. 

C 

NDIFF = IA8SU-K) 

r 

Vn 

C SCHEME A3 ANDONED IF DIFFERENCE GREATER THAN BANDWIDTH OF 

C PRIVIOUS SCHEME, NEW SCHEME STARTEO. ‘ 

C • ‘ 

IFIND IFF.GE .MINMAX) GO TO 60 
IFINDIFF.GT .MAX) MAX =NDIFF 
40 CONTINUE 

I F { K. EQ. NODES )G0 TO 50 ‘ 

45 1=1+1 
’ GO TO 30 
50 M INMAX = MAa 
DO 55 J = l, NODES 
55- JNT ( J } = JOINT! J) 

60 CONTINUE 
RETURN 
ENO 



000002 

000002 

000002 

000002 

000002 

000004 

000013 

000016 

000017 

000022 

000024 

000026 

000027 

000031 

000034 

000036 

000040 


000042 

000055' 

000055 

000077 

000122 

000122 

000123 
000123 
0001 25 
000130 
000153 

000200 

000200 

000201 

000201 

000.204 


SUBROUTINE SAPOUT 

COMMON LMENTS, JT(4Q00) , ME MJT { 8000 ) , JM EM (1000) » JNTC 1000 1 

COMMQN/CDNTR/ NELG, I TYP E ( 5 I ,N EL ( 5 ) , NODES , IPUNCH 

COMMON/ JUNK/A (1000, 61 

COMMON/UN I T / I N , I T f IP 

DO 10 1= 1, NODES 

WRITE (IT, 12) I , JNT{ I ) 

10 CONTINUE 
LMENTS = 0 
WRITE(ET,14) 

DO 20 I 1 = 1 , NELG - 

N = ITYPEtm 

NE= N EL (I II 

DO 21 JJ =1 , NE 

I = JTI J J+LMENT o , 

J = JT( JJ+LMENTS +1000) 

N 1= JNT(I) 

NJ= JNT(J) 

C 

C OUTPUT NODE CONNECTIONS ACCORDING TO TYPE. 

C 

GO TO ( 1,2, 3, 3., 5,3,7), N 
C 

1 CONTINUE 

WRITE! IT, 102) JJ,NI,NJ, (A(JJ,L) »L=1 ,2) 

IF(.IPUNCH.GT.O) 

*W RITE l IP, 10 2) JJ,NI»NJ, (A(JJ,L) ,L=1,2) 

102 FOR MAT! 31 5,2410) 

GO TO- 21. . 

C 

- 2 CONTINUE ■ 

K, = A ( J J, 1 ) 

NK=- JNT C<) 

WRITE (IT, 104) JJ,NI,NJ,NK, ( A( JJ ,U ,L=2» 61 
IF( IPUNCH. GT.O) - 

*WRITE ( IP, 104) JJ,NI,NJ,NK, ( A( JJ ,L ) ,L=2, 6 ) 

104 FORMAT! 41 5 , 5A 10) 

GO TO 21- -- 

C 

3 CONTINUE 

K = JT( JJ+LMENTS +2000) 

L = JT( JJ+LMENTS +3000) 



000206 



IF(L.EQ.O) L=K 

000210 



NK= JNT(K) 

000212 



NL= JNTIL) 

000214 



WRITE! IT, 106)JJ,NI,NJ,NK,NL»< A( JJ,L),L=L,5> 

000242 


106 

£= 0R M AT[ 51 5, 5410) 

000242 



I F { IPUNCH.GT.O) 




*WRITE (1°, 106) JJ,NI,NJ,NK,NL,( AI JJ,L) ,L=l,5) 

000271 

r ' 


GO TO 21 

000272 


5 

CONTINUE 

000272 • 

r 


GO TO. 21 

0.002 73 


7 

CONTINUE 

000273 . 



WRITE(rT,108)NI,NJ,(A( JJ,L).»L=lt6) 

000313 



" I F ( IPUNCH.GT.O 1 




-WRITEUP, 108)NI ,NJ,(Al JJ,L) ,L=1,6) 

000334 


108 

FORM AT ( 21.5 , 6A 10 J 

000334 


21 

CONTI NUE 

000337 



LMSNTS = LM ENTS +NE 

000340 


20 

CONTINUE . - 

000342 


12 

FORMAT ( 16 x", 15 , 13X, 151 

000342 


14 

FORMAT! 1H1 , //,4X, 17HNEW ELEMENT CAROS, // ) 

0003 42 


.30 

FORMAT ( 515) 

000342 



RETURN 

000343 



END 



APPENDIX C 
SAMPLE BANSAP OUTPUT 
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CQ CO CO CO CO 


BBBBB . 
B 
B 

BBS 

B 

B • 

BBBBB 


AAAAA 

N 

N 

SSS5S 

AAAAA 

A A 

NN 

N 

s 

A A 

A - A 

N 

N N 

s 

A A 

AA AAAAA 

M 

N N 

ssss 

AAAAAAA 

A A • 

N 

N N 

s 

A A 

A A 

N 

NN 

s 

A A 

A A . 

N 

' N 

sssss 

A • A 


PPPPP 

P P 

P P 

PPPPP 

P 

P 

P 
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INPUT DATA 


NUMBER OF 

ELEMENT TYPES = 

1 


NUMBER OF 

NODAL 

POINTS = 

9 


PUNCHED ELEMENT 

CARDS 

0 



.EQ. 

0 NO 




.EQ. 

1 YES 



ENT GROUP 1 

HAS 

8 'ELEMENTS 

OF TYPE 

3 

ELEMENT ' 

- I . 

J 

K 

L 

1 

1 

9 

8 

-0 

2 

i 

2 

g 

-0 

3 

2 

4 

9 

-0 

A 

2 

3 

4 

-0 

5 

8 

6' 

7 

-0 

6 

8 

9. 

6 

-0 

7 

9 

5 

6' 

-0 

8 

9 

4 

5 

-0 
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NODE J MEM MEHJT 


2 

4 3 

3 5 - 

4 

9 5 

6 7 

2 4 6 5 


ORIGINAL 'BANDWIDTH = -9 


NEW BAN DWIDTH= 4 
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OLD NODE NUMBER NEW NODE NUMBER 


1 

2 

3 

4 

5 

6 
7 
8 . 
9 


4 
2 
1 
3 
6 
8 
9 
7 

5 
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NEW ELEMENT CAROS 


1 4 5 7 7 

2 4 2 5 5 

3 2 3 5 5 

4 2 13 3 

5 7 8 9 9 

6 7 5 8 8 

7 5 6 8 8 

8 5 3 6 6 
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Table 1. Example node connections determined 
in subroutine SETUP. 


' Node 

Number of 
Connected 
Nodes 

Connected 

Nodes 

1 

3 

9, 8, 

2 


2 

4 

1, 9, 

4, 

3 

3 

2 

2, 4 



4 

4 

2, 9 / 

3, 

5 

5 

3 

9, 6, 

4 


6 

4 

8, .7, 

9, 

5. . 

7 

2 

8, 6 



8 

4 

1, 9, 

6, 

7 

9 

6 

> 

CO 

2, 

4/ 6, p 
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NEY/ NODE NUMBERS WROTE RELATED NODES ARE ASSIGNED NEW NUMBERS. 



Table 2i Trail Numbering Schemes used in OPTNUM 


OLD NODE NUMBER AT WHICH ORIGIN OF NEW NUMBERING SCHEME IS SET. 


! X | X I 


2 — i - — 5 
I \ I \ t 
6 — 3 4 

1 — 7— B 

DIFF- S 

ABANDON 

SCHEME 


I IX I 


DIFF -2' 

4 r — 2 — 1 

j X f X I 

— 5 — 3 

IX |X I 

DIFF -'3 
4 — j — 1 

•liU 

'Hit 


4 — 2 — — 1 

I \ I \ I 
7 — 5 — 3 

II '1J 


4 — 2 — 1 
l\l\i 

7—6 — 3 

'1J1J 



DIFF » 4 

ABANDON 

■SCHEME 



|M1 

Hill 

DIFF -3 
ABANDON 
SCHEME 



IX |\J 

1151! 


DIFF -4 
ABANDON 
SCHEME 


DIFF -2 


i\l\| 

I — 3 v — ' 

Hill 

DIFF -4 


nil 
1111' 
DIFF -3 
ABANDON 
SCHEME 


— ABANDON 


SCHEME 


OIFF- LARGEST DIFFERENCE j 
BETWEEN ANY TWO RELATED NODES, 
































Table 3. Summary. of applications of BANSAP 




Number 

Number 

Old 

New 


Element 

of 

of 

Band- 

Band- 

Structure 

. Type 

Nodes 

Elements 

width 

width 

Sample problem • 
Figure 1 

Membrane 

9 

8 

. 9 

4' ■ 

Wagonwheel 
Figure 6 

Truss 

9 

16 

9 

6 

Ship tower 
Figure 7 

Beam < 

25 

65 

12 

9 

Shear panel specimen 
Figure 8 

Membrane 

59.5 

554 

406 

35 

Bolted joint specimen 
Figure 9 

Membrane 

398 

349 

168 

28 












































START • 



STOP 


Figure 5. Flowchart of SAP IV 
BANSAP Preprocessing 
Program. 







(a) Original Scheme 



(b) Renumbered Scheme 


Figure 6. Wagonwheel Truss. 
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Figure 8. 


Shear Panel Specimen. 
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TRIANGULAR AND QUADRILATERAL MEMBRANE ELEMENTS 



Figure 9; Finite Element Mesh for Composite Bolted Joint Specimen. 


FEMESH: 


A FINITE ELEMENT MESH GENERATION PROGRAM 
BASED ON ISOPARAMETRIC -ZONES ' 


By 

Zoa C. Lane' 
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FEMESH: A FINITE ELEMENT MESH GENERATION PROGRAM 

BASED ON ISOPARAMETRIC ZONES 

By 

Zoa C. Lane 
INTRODUCTION 

Finite element analysis programs greatly facilitate the 
determination of deformations and stresses in structures. A 
major inconvenience in utilizing this analysis technique is the 
large amount of input data required by the computer programs. 

This data includes,, in addition to material characteristics, 
the' node .numbers defining the elements and the spatial, coor- 
dinates for each node. 

Current mesh generation methods inciude for- simple problems 
data preparation - by hand-, and for more complex problems, the coding 
and .executing of FORTRAN mesh generation programs which generate, 
data for a- general structure. 

W._R. Buell and W.A: Bush surveyed some techniques, used in 
current mesh generation schemes (ref. 1). The techniques pre- 
sented by Buell and Bush are: a straight line interpolation tech- 

nique, a sides and parts technique for axisymmetric structures 
electro-mechanical techniques for two- and three-dimensional 
structures, and a simplified finite difference technique and 
equipotential technique for- general structure shapes. 

The advantages of general structure mesh generation programs 
(ref. 1) are: (1) reduced cost due to reduction of man hours 

and computer time needed to generate and check data; (2) reduced 
number of errors; (3) insured regularity of finite elements; and 
(4) application to a variety of structural shapes. 

O.C. Zienkiewicz (ref. 2) utilizes a technique involving 
the mapping of isoparametric quadrilaterals from a natural to 
a cartesian coordinate system in an automatic mesh generation 
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scheme for plane and curved surfaces. This scheme is applicable 
to non-quadrilateral structures if the structure is divided into 
quadrilateral regions. Zienkiewicz * s technique for mesh genera- 
tion was used by S.J. Womack (ref.- 3) as a preprocessor for TEXGAP, 
a finite element program for the analysis of two-dimensional 
linearly elastic plane or axisymmetric bodies (ref. -4) . 

The objective of this study is - to utilize the technique 
developed by Zienkiewicz in a mesh generation scheme for two- 
dimensional planar surfaces. Presented in this paper are a 
description of the mapping technique, a description of the 
computer program; and three examples of meshes' generated by 
the program. A set of user instructions and a listing of the 
program are included in the appendices! 

■ INTERPOLATION FUNCTION TECHNIQUE' FOR FINITE 
ELEMENT GENERATION. 


The algorithm used by Zienkiewicz to ' map an isoparametric 
quadrilateral is the displacement' interpolation equations used 
in isoparametric, finite elements (ref. 5) . The interpolation,' 
equations for quadratic bounded surfaces (which, are listed in 
table 1), are a function of a set of dimensionless .coordinates, 

C and r, which define a natural' coordinate system. 

In the natural coordinate system (fig. 1) , a planar* surface 
is represented as a square whose dimensions are 2x2 units and 
whose center is at the origin. To map a surface into the ■ 
cartesian coordinate system, eight boundary points- (x^ and y^) 
and the £ and r values of each grid point on the surface 
to be mapped are substituted in the displacement -interpolation 
functions; the resulting values are the cartesian coordinates 
of the grid points. 

A mesh is generated by dividing, the square into the desired 
number of subdivisions , calculating the £ and rj coordinates 
for each grid point, and mapping each point to the cartesian 
coordinate system. A graduation of a generated mesh is obtained 
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by offsetting the midside node from the midpoint of a side of- the 
quadrilateral (fig. 2) . The generated elements will vary' in size 
along that side; smaller elements will be in the direction of 
the offset. - 

Meshes for complex structures are generated by dividing the 
structure into quadrilateral zones. The mesh for each zone is 
generated independent of other zones. Connection of zones is 
accomplished by eliminating node numbers and coordinates which 
were duplicated on zone boundaries. 

PROGRAM FEMESH 

Program FEMESH is a FORTRAN IV code for generating finite 
element data for two-dimensional planar surfaces. The algorithm 
used to generate the node coordinates is based on the displace- - 
ment interpolation functions (table 1) described in the preceding 
paragraph. 

Input data for FEMESH includes a title, the number of zones, 
.the total number of zone, nodes, the number. of zone node coordinates 
to be read- from cards’, the first node and element numbers, a - list 
-Of- the _ eight nodes which define a zone, the dimensions of the 
desired mesh of each zone, and the zone node coordinates. 

A zone is a quadrilateral region whose geometry is defined 
by eight -zone nodes. (Zone nodes are used only in the input . 
definition- of the geometry; they are. not included in the generated 
mesh.) The zone nodes are listed in counter-clockwise order. 

As indicated in figure’ 3, the first node identifies a corner of 
the quadrilateral. The second, fourth,’ sixth,- and eighth nodes 
are referred to . as midside nodes. If a midside node does not 
lie on the midpoint of a side, a graduation of the mesh results. 

The general flow for the mesh generation program, FEMESH, 
is shown in figure 4. As indicated, the mesh for each zone 
is generated separately.- The first step in the mesh gener- 
ation scheme is to determine if the coordinates of the midside ‘ 
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nodes are defined (i.e., if their coordinates were supplied 
by the user) . If the coordinates are not defined, the midside 
node is assumed to lie at the midpoint of a linear line seg- 
ment. The second step is to determine if either of the four 
sides of the zone is connected to a zone for which a mesh 
was previously generated. If a side is connected to such a 
zone the node numbers and the x and y coordinates which 
have already been generated are used. The remainder of the 
mesh xs then generated; This process is repeated until .the 
meshes for all zones are generated. 

The output of program FEMESH includes a listing of the 
elements, their four node numbers and the node coordinates. A 
plot of the mesh is also generated. 


APPLICATIONS 

Three finite element mesh , generated by 'FEMESH are presented 
in this section. The first example is a sample problem illus- 
trating the input and output of program FEMESH. ' The second is 
a quarter section of. a shear panel. • The third is a half section 
of a bolted joint specimen. 

. . The first example is a simple structure originally used to 
validate the ability of FEMESH -to properly connect zones. The 
structure (illustrated in fig. 5a) is divided into three zones; 
The- eighteen zone nodes .are labeled arbitrarily and illustrated 
in figure -5b. Figures 5a and 5b represent the input required,- 
by program FEMESH to. generate the mesh' illustrated -in -figures 
5c and 5d. Figure 5c illustrates. the node numbers, and 5d 
iiinstrates the' element numbers. 

The input data for this problem is tabulated in table 2 ' 
(see Appendix A for user' instructions) . The data includes a 
title card, a control card, three zone description cards, and 
eight node coordinate cards. The control card specifies the 
number zones (3)., the number of zone nodes (18), the number 
of zone node coordinate cards to be read (8), the first node 



number (100) , and the first element number (1000) . A typical- 
zone description card lists the eight zone nodes defining each 
zone and the size of the finite element mesh to be generated. 

The tabulated output for this problem appears in table 3.‘ 

The output includes the input data, the element number, the 
four node numbers which define each element, and the cartesian 
coordinates of each node. 

The shear panel illustrated in figure 6 is divided into four 
zones. The zones were established in such a way that the straight 
and curved segments of the corner fillets are assigned to different 
zones in order to obtain a closer approximation of the true 
boundary shape. 

The mesh dimension for zone .1 is 20 x 20, for zone II is 
20 x 3, for zone III is 20 x 20, and zone IV is 3 x 30.. To 
avoid the generation of long, narrow rectangular elements, the 
midside nodes 2, 8, 9, and 14, 15, 16 are moved away from -the 
midpoint of the line segment toward the fillets. The input data 
is summarized in -table 4. The' output is -illustrated in figure 7. 
Because of the large number of generated elements, the output is ' 
not listed in tabular form; it is represented graphically by a' - 
'computer plot of the generated" mesh. The generated mesh is com-, 
posed of 5-74 nodes-- and 52-0 elements. 

The mesh for one-half ' of a bolted joint specimen was 
generated by. dividing the specimen into 15 zones as illustrated 
in figure 8. The input data for this problem (table 5)- consisted 
of 58 data cards, including 15 zone description cards, and 41 node 
coordinate cards. A graduation of the mesh of zones II, III,. IV, 

V, VI, VIII, and IX was used to obtain a uniformity in the shape 
of the generated elements. The generated mesh, which is illus- 
trated in figure 9, consists of 378 elements and .435 nodes. 

CONCLUDING REMARKS 

Program FEMESH, a FORTRAN IV code, has been developed to 
generate a finite element mesh for two-dimensional, planar , 
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surfaces. The algorithm used is the displacement interpolation 
functions which were developed for mesh generation by Zienkiewicz, 

A structure may be subdivided into a maximum of 15 zones. 

The maximum mesh, for each zone is 2'4 x 24 elements (or 25 x 25 
node points). FEMESH will compute a maximum of 4000 node points, 
and output the node numbers and their coordinates and the element 
numbers and their four identifying node numbers. A simple plot 
of the finite element mesh is also generated. 

Presented in this paper is a description of the technique 
used in the mesh generation scheme, a description of program 
FEMESH and examples of the mesh generated for three problems. 

User instructions and a listing of the program are included 
in the appendices. 



APPENDIX A 

USER INSTRUCTIONS FOR FEMESH 
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USER INSTRUCTIONS FOR FEMESH 


Program FEMESH generates isoparametric finite element meshes 
for two-dimensional planar surfaces . The input required by the 
program consists of four types of data cards:- a title card, a 
control card, zone description cards, and node coordinate cards 
(fig. Al) . 


TITLE CARD (Format 10A4) . 


Column 

Variable 

Description 

1-40 

TITLE 

Heading for output 

, CARD 

(Format 615) : 


Column 

Variable 

Description 

1-5 

IZ 

Number of zones (IZ < 15) 

5-10 

NT 

Total number of zone node^ 

11-15 

- NI- 

Number of zone node coordinates to 


- - 

be read as input on cards 

16-20 

.INODE 

First node number to be assigned to 



generated mesh 

21-25 

IELM • 

First element number to be assigned 



to generated mesh 

26-30 . 

. . IP 

Punch indicator: 0 will not punch ■ 


1 punch 


A zone. is a quadrilateral with either linear or curved line 
segments. The geometry of the zone is defined by 8 zone nodes 
.whose coordinates are supplied by the user (see node coordinate 
card) . 

The values of NI and NT may differ due to the ability 
of the program to linearly interpolate to define the coordinates 
of the midside node if those coordinates are riot supplied by the 
user. Midside nodes are those zone nodes which lie between two 
corner nodes. It is not necessary that a midside node lie- at 
the midpoint of a line segment. 
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ZONE DESCRIPTION CARD (Format 


Column 

.'Variable 

1-5 

NODE 

(1,1) 

6-10 

NODE 

(1,2) 

11-15 

NODE 

(1,3) 

16-20 

NODE 

(1,4) 

21-25 

NODE 

(1,5) 

26-30 

NODE 

(1,6) 

31-35 

NODE 

(1,7) 

36-40 

NODE 

(1,8) 

41-45 

M 


46-50 

N 



1015) : 

Description 

Zone nodes defining zone 
geometry 

I is the zone number 


Number of subdivisions along the 
side defined by 1st, 2nd, and 
3rd zone nodes 

Number of- subdivisions along the. 
side defined by 3rd, 4th, and 
5th zone nodes. 


Zone numbers are determined by- the order of the zone des- 
cription cards. The- first zone description card is assigned the 
number one, the second is assigned the number two, etc. 

The interconnectivity of zones is indicated by assigning 
a negative magnitude to zone nodes which lie on a side connected 
to a zone with a smaller zone number. For example, if 4 zones 
are connected as shown in figure A2,.then the first eight values 
of the zone description cards should be: 


Card 

1: 

1 

2 

3 

7 

11 

.10 

9 

6 

Card 

2: 

-3 

4 

5 

8 

13 

12 

-11 

-7 

Card 

3: 

-11 

-12 

-13 

16 

21 

20 

19 

15 

Card 

4: 

-9 

-10 

-11 

-15 

-19 

18 

17 

14 


A side which is divided into M subdivisions must not be 
connected to a side divided into N subdivisions unless the 
values M and N are equal. 
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NODE COORDINATE CARD (Format 15, 2F10.5): 


Column 


Variable 


1-5 

5-15 

16-25 


Node number 
x coordinate 
y coordinate 


This card may be omitted for any midside node which lies on 
a straight line if a graduation of the mesh is not desired. 

A graduation in the mesh occurs .when the midside node is 
offset from the midpoint of the line segment. The smaller 
elements will be in the same direction as the offset. 

Due to a restriction in the FORTRAN coding, a midside node 
should not be assigned the coordinates (0,0) if the line seg- 
ment is 'not a straight line. 
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ONTINU ATION 


FORTRAN STATEMENT 


T I 9 10 11 12 ,11 14 |S l» 12 II J? Jl 22 21 Zi 11 2A 21 28 22. 20 11 12 1J 14.11 It 12-11 11.40 11-12 _0.44_4i_4*_4Z_j 


iLiXiTiLEt. i .j — i. -i- 


.J L. i 1 1 — 


. X J J 1~ 1- 


ggg||£j 1(93358 5 

^5l— HMM 


£ 0.N,T>Ka.L^C t AK& 


XMOiOE .. iXjEit-M -r-i-_.iX-iP.i_. _.i — i— i _i..j — i — i — i — i_i — 

■ I.' ■ ■ ■ ■ - • ■ F'iCJiKAlATi 


,j. i. j i. t i„ i r i_..i 1. 1. •!. x — i . . j i .i_ j -i~- • — i_.j . j — J. — i — 


I.. 1. J. i...l -X- L 1 L-d _L. L..L .1. . I . J. X - L-J. -1— J — L_ I ! I I I 1 ^ 1 

1 Z0.MF DfSC^iPTIfl 


el ,Ni£i.DiE 


L _X,..4 .X „ X .1 


. J L-. I L, 


Ni a <£M .RxtuX 


I 'MATE 1 c ,0, < y,R,b|X l NAT,g| ^0^MAT 1 ^y^iZJFa^J^ 

II I I L_ 1_J I l__L. .J 1 1 ' „1 — J — J 1 — ,1 -- J — J-' — I * 1 1— 

1 i \ . * r 

J I » J .J l X..I X. X L.J L.,.I._J r I.J_l...X..l .1 .1—1 ‘ I— I J 1... 1 ...I ... 


L-. 1 L I I -J . l—.i — -X._J. 

I L. J I I 1 L.: U..I L—J — I 1 — 

.L J L J t — L_J 1 ' J — !_ 1. 1 1-X.X_ 


. Figure ,A1. 1 Input Data Formats for, FEMESH. 


ui 











Figure A2. Simple Structure to Illustrate 
Zone Node Input Data. . 
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APPENDIX B 


FORTRAN LISTING OF MESH GENERATION PROGRAM, 
(LRC, CDC-6600 COMPUTER VERSION) 


FEMESH 



000003 

000003 

000003 

000003 


000003 

000004 


000005 

000006 


000007 

000007 

000007 


PROGRAM FEMESH! INPUT, OUTPUT, TAPES* INPUT, TAPE 6=0UTPUT, PUNCH) 

C 

C 

c . , 

C PROGRAM FEMESH CODED BY Z. C. LANE MAY 31, 1975 

C 

C 

C '.PROGRAM FEMESH GENERATES FINITE ELEMENT DATA FOR TWO DIMENSIONAL 
C PLANAR SURFACES. STRUCTURES MAY BE SUBDIVIDED INTO AS MANY AS 15 

C UUADRI LATERAL ZONES* THE MAXIMUM MESH DIMENSION IS 24 X 24 SUBDIVISIONS. 

C 
C 
C 

C * NEGATIVE ZONE NODE NUMBERS FOR A ZONE IDENTIFIED BY A NUMBER *j*N** 

C 1DNICATES’ THAT .THE NEGATIVE NODE -IS CONNECTED TO A ZONE WHICH IS 

C IDENTIFIED BY A NUMBER LESS THAN **N**. 

C 

C , • ■ 


** 

INODE 

FIRST NODE NUMBER 

' ## 

M,N 

NUMBER 

OF 

SUBDIVISIONS 

** 

I CTZ 

CURRENT ZONE NUMBER 


NT 

TOTAL 

NUMBER OF 

INPUT NODES 


I Z 

TOTAL NUMBER OF ZONES 


N I 

■ NUMBER 

OF 

NODES 

TO BE READ 

** 

IP 

PUNCH WHEN IP=T 

, <c* 







c 

c 

C 

DIMENSION A( 8), TITLE ( 101 .XNCDEIB) ,YN00£(8, 
DIMENSION X1(78),YU78) 

COMMON ICTZ»IZCNE(15»10J ,N00£ (15 , 25 ,25 ) , TEMPI25 ) 
COMMON X2(4002i,Y2(4002) 

C 

c 

CALL PSEUDO 
CALL LEROY . 

C 

C * DETERMINE, INPUT - OUTPUT DEVICES 
IN = 5 
IOUT = 6 
C 

c 

1 FORMAT (615) 

2 FORMAT ( 101 5) 

3 FUKMAT(I5,2F10«5) 



CC0G06 

000006 

000006 

000006 

CCC0C6 

000006 

000006 

000006 

( 


0000 06 
000010 
OCOOll 
000012 


000014 

000021 

000041 


000060 

0C0066 

000074 

OCOIOO’ 

0.00121 

000125 

000127 

000140 

000152 


000155 

000156 

000157 


3 F0RMAT(I5»2F10.5) 

4 F0RMAH1H ,12HNQ. OF ZONES, I 2,//) 

5 FORMAT ( IH » I 3 , 1 X ,8 15 , IX , 2 14 ) 

■6 .FORMAT ( 1H ,I3,2X,F7-3,2X,F7.3) , \ 

7 FORMAT { 1 0A4) 

8 FORMAT! 1H 1 , l 0A4 /// ) 

9 FORMAT ( 1 H , 4HZ0N E, 15X, L OHZON E NODES, 1 8X , IHM , 4X , 1HN , /2X, 3HNQ. ,4X»IH 
1 1,4X, ; 1H2,4X, 1H3,4X,1H4,4X,1H5,4X,IH6 ,4X , 1H7 ,4X , 1H8/ ) 

12 , FORM AT !// / , 1 X, 4 HNQDE , 4X, 1HX,'8X, 1HY ) 

C ■ 

C 

C 

C 1 INITIALIZE XI AND . Y l TO BE FILLED FROM DATA READ OFF CARDS 
DO 34 1=1,78 
Xl( I) = 0. 

■ V 1 ( 1 I =0. 

34 CONTINUE 
C 

c 

c 

C * READ INPUT 

■ READ! IN, 7 ) Yl TL6 

READ! IN, 1 HZ, NT, NI, INO.DE, I ELM, IP . 

READ! IN, 2 I! ( IZQNEU, J), J=1 ,10 J ,1 = 1, IZ) 

C 

C * WRITE INPUT 

WRITE! 10UT ,8 >Tl TLE 
WRITE { IOUT ,4 I IZ 
WRITE! IOUT ,9} 

WRITE (IOUT ,5 I ( I, ( I ZONE! I, J), J=1,10|, 1=1, IZ » 

WR ITE! IOUT ,12) - 

c' 

CO 10 J= 1 , Nt 

READ! IN, 3) I ,X1(I>,Y1{J) 

WRITE! IOUT ,6 ) ItXHIJtYim 
10 CONTINUE. 

c 

C * SET COUNTER OF NODE NUMBERS, ICTN AND Z0NE,ICT2 
ICTZ = 0 
IC'TN = INODE 
NCOR = 0 ‘ 

C 



000160 
000161 
000162 
000163 
GOO 1 6 5 
000166 
000167 
000170 
000175 
0C0177 
00020 


00020 3 
C002C3 


000205 
0C0206 
0002C7 
0002 10 


000211 

000212 

000216 

000220 

0C0221 

000223 

000225 

000226 

000230 


000231 

00023T 


00 35 1 = 1,1661 

*21 I) = 0. ' 

Y2U ) = 0. 

35 CONTINUE 

CC 1010 Ifl.13 ■ 

CO 1009 J=l, 25 
DO 1008 K=1 ,25 
NODE ! I, J,K> = 0 

1008 CONTINUE 

1009 CONTINUE 

1010 CONTINUE 
C 

G 

1000 CONTINUE 

ICTZ = ICTZ + l 
C 1 
C 

C , SET I DNI CATC1RS TO ZERO 
C 

151 = 0 

152 =:0 

1 S3 =0 
IS4 = 0 

C 

C , * PULL FROM THE IZONE ARRAY THE ZONE NUMBERS AND THE ZONE COORDINATES 

C ‘ 

DO 20 1=1,8 . 

IC = IABSI IZONE 1 ICTZ, ID 
XNODE(I) = XIII Cl 
YNODEU) ,* YlltCl 
20 CONTINUE 

K •= IZONE! ICTZ, 9 1 
N = IZONE! ICTZ, 101 
'NN f N +1 
MM = M + 1 , 

C * 

C’ 

C 

C TITLE' THE OUTPUT FCR THIS ZONE 

C 

WRITE! IQUT.lllICTZ 

11 FORMAT (1H1 ,28HCALCULATIONS FOR ZONE NUMBER , 1 4/// 1 


oo 

o 




C 

• c 

IF' NO VALUE IS GIVEN FOR MIDPOINTS* ASSUME A STRAIGHT 


c 

r 

TFE MIDPOINT'.’ ‘ 

000237 


CO 30 1 = 2*8, 2 

0002 A 1 


I F ( XNODEI 1)130,25, 30 

0002 A 2 

25 

IFIYNODEC I ))30, 26,30 

000244 

26 

K * 8-1 

0C0246 


IF(K)30,27,28 

000247 

28 

XNODEI I) = ( XNODE (I + D+ XNODE (1-11 )/2. 

000253 


VNODE (I) = ( YNODEU + l)+YNODE(I-l) )/2. 

000256 


GO TO 30 

000256 

27 

XNODE 1 1 ) = ( XN0DEU)+XNaDEI7})/2- 

000262 


YNODEII) = ( YNODE(l) ♦YN0DEI7 ) )/2.' 

000264 

30 

CONTINUE 

000266 


WRITE I I OUT ,31 ) I I , XNC-CEI N ,YNOCE{ 11,1 = 1,8) 

000303 

31 

FORMAT ( 1H ,10HZ0NE NODE , 14 , F 15. 2 , F25. 2) 

000303 


WR I TE ( I OUT ,32) 

000307 

32 

C 

C 

FORMAT <//* IX, 8HN0DE. NO. , 5X, 1HX,8X , 1HY/ ) 


C 
C 
' C 

* IF ZONE 'NUMBER IS ONE, FILL NODE ARRAY AND SKIP TO X, 


C 

r 

, CALCULATIONS’. 

000307 


IF! ICTZ-1)199,1S0,2G1 

CC03I2 

190 

CO 192 J=i ,NN 

000314 


DO 191 1 = 1, MM 

0GO315 


NODE (1,1 ,J) =• ICTN 

000322 


ICTN = ICTN * 1 

000323 

151 

CONTINUE 

000325 

192 

CONTINUE 

0CO330 


GO TO 3000 

000330 

1.99 

WRITE! IOUT ,200) ICTZ 

000336 

200 

FORMAT! 1H1 .28HERR0R ... ZONE NUMBER ICTZ =.141 

000336 

201 

CONTINUE 


DETERMINE WHICH SIDES ARE CONNECTED 
FILL THE NODE ARRAY 


LINE AND CALCULATE 


Y COORDINATE 



000336' 

000 340 

000343 

000345 

000353 

CC0355 


000356 

000356 

000360 

000363 

000365 

000375 

000377 


000400 

0C04C0 

OC04G2 

000405 

Q004Q7 

0C04L7 

000421 


000422 

000422 

000424 

CCC427 

000431 

000437 

000441 

0C0442 


C 

C. SIDE ONE 
C 

. IF! I ZONE l ICTZ,2.)'>212,220,220 
212 , CALL F INDI 1 , 3 , 1 1 
CO 213 1=1 ,MM 
' NODE I I CT Z, 1,1) *■ TEMP ( I ) 

213 CONTINUE- 
I SI = 1 

c 

c 

C . SIDE 2 

C 

220 CONTINUt 

. IFIIZCNEI ICTZ, 41 )222 ,230,230 
222 CALL FIND (3, 5, 2) 

CO 223 I =,1 ,NN 

' NODE ( ICTZj MM, I ) = TEMPI I ) 

1 223 CONTINUE 
IS2 = 1 
C 
C 

C SIDE 3 
C 

230 CONTINUE 

IF! I'Z ONE I I CTZ, 6) 1232,240,240 
. 232 CALL F I ND I 5, 7, 3 1 
\ CO 233 1=1, MM 

NODE! ICTZ,I,NN) .= TEMPII 1 
• 233 ' CCNTINUF ? 

■ ■ IS 3 = 1 

C . 

c 

C SIDE 4 

C 

240 CONTINUE 

I F ( I ZONE 1 1 CTZ, a 11242,250,250 

242 CALL FIND! 7, 1 ,4) 

• CC 243 1=1, NN 

NODE { ICTZ , l'» I 1 = TEMP ( M 

243 CONTINUE 

■ IS'4 = 1 1 
250- CONTINUE 



000442 

OC0444 

000445 

000453 

000453 

000461 

000462 

0C046 3 
000463 
000467 
CCC467 
000472 
0CC474 


000474 

000476 

000477 

C005CI 


. 000503 
CC05C4 

0C0505 
000506 
0C051C 
000513 
0005 1 5 
000517 


C 

C 

C 

C FILL NODE ARRAY - JUMP THOSE PQS I T IONS. A LL READY FILLED 
C 

C > 

DO 320 J = 1-,NN 
CO 310 1=1 , MR 

I F ( NODE ( IC TZ , 1 , J) ) 315,300,310 
300 CONTINUE 

NODE I ICTZ , I , J 1 =' IC Tl\ 

I CTN = ICTN' + 1 ' 1 
GC TO 310 
C . 'ERROR 

315 CONTINUE 
WRITEtiaur ,316) / 

316 FORMAT ( 46H NODE NO FOUND IN ST. N/ . 300-320 LESS THAN ZERO) 

310 CONTINUE 

.320 CONTINUE 
3000 CONTINUE 
C 
C 

C * COMPUTE THE X-Y COORDINATES, OM I TT ,PRE VI QUSLY COMPUTED SITES. 
C 

C 4, DC, AND CN ARE THE INCREMENTAL ' VALUES IN THE M AND N 

C DIRECTIONS RESPECTIVELY. 

, PM = M 
PN = N 
' CC = 2 ./ RM 
ON = 2. /RN 

C • - kRITEI I OUT ,33) DC, ON 

C 33 FORMA T 1 1H ,4HDC = . ,F5 . 2 ,4X ,4HDN= ,F5.2) ' 

C 

CCC = -I.- 
CCN = -1 . 

C 

DO 810 J=1 ,NN • 

IF! J— 1)731,730, 731 

730 IF! I S 1- 1 1 733 , 51 ,733 

731 IFU-NNJ733, 732,733 

732 IF! IS 3- 1)733,5 1,733 

733 CONTINUE 
C 


co 

u> 



000517 

.c' 1 

CO 800 1=1 , MM 

000521 


I Ft I-MM >741, 734',- 741 

0C0523 

734 

r 

I F l IS2-1 J739‘,50,739 

0Q0526 

741 

I F t 1-1)739, 742, 7 39' 

0CC53G 

• 742 

1 F ( IS 4- 1)739, 50, 739 

000532 

739 

CONTINUE 

000532 


SI = l.-CCC 

GCQ534 


S2 = l.-CCN ' 

000536 


S3 = CCC+CCN-l. 

000540 


S4 = l.+CCC 

000541 


15 = l.+CCN 

C00543 


Ml) = l./4.*Sl*S2# ! -CCC-CCK— 1 . ) 

000552 


M2) = l./ 2. *31*52*54 

0C0554 


M3) = l./4.*S2*S4*(CCC-CCN-l.) 

000563 


M4 ) = l./2.*S4*S2*S5 

0C0566 


M 5 ) = ■ 1. /4.*S3*S4*SS 

000571 


A<6) = l./2.*Sl*S4*S5 

000575 


A ( 7 ) = l./4-.*Sl*S5*(-CCC«-CCN-l.) 

CC06C4 


A ( 8 ) = l./2.*Sl*S2*S5 

000610 


NCOR = NCOR + 1 

000611 


CO 45 K=1 ,8. 

000613 


X2INC0R) = AIK) *,XNODE!K), + X2INC0R) 

0C0617 


V21NC0R) = A ( K) * YNCDE t K ); + Y2CNCCR) 

000622 

45 

CONTINUE • ’ . 

000624 


WR I T E ( IOUT ,2000 ) NC GR , X2 ( NCOR) ,Y2(NC0RJ 

OC0635 

2000 

P0RMATU5,2X,F10.5,2X,F10.5) ■ 

000635 

50 

CONTINUE 

0C0635 


CCC = CCC+DC 

0C0637 

800 

CONTINUE' 

000642 

.51 

CONTINUE , 

0'CC6 42 


CCN = CCN + OK, 

000&44 


CCC = -1. 

0C0646 

810 

CONTINUE 

000650 


WRITEUOUT ,54) 

000654 

54 

FORMAT!///, 1H , 1 5HN0DE NO. MATRIX/) 

CCG654 


CO 53 J=1,NN ! 

0006S6 


WRITE! 10 UT ,52 M NODE IICTZ,I»J) ,1=1, MM), 

CCC6 73 

52 

FORMAT ! IX , 26 15) 

000673 

53 

CONTINUE 

000676 


IF! ICTZ-IZ)1000,4000C, 40000 

CC07C0 

40000 1 

CONTINUE 


00 



CC07GC 
0 Q 0 7 0 4 
CC07CA 
000706 
0C0710 
000711 
000713 
0C0714 
000716 
0C0717 
CCC724 
000731 
000735 
000742 


000757 

000777 

000777 

CC0777 

001001 

OCIOC3 

001006 

001010 • 


001013 
001017 
CCl 017 
001021 
0C1032 
CCl 032 
001034 


C ... 

C * LIST THE ELEMENT NUMSERS AND OEFINING NODE NUMBERS . 

C ■ 

C 

' WRITE ( IOUT ,900 1 ' 

. 9D0 'FORMATt 1H1,7HELEMENT ,8X,1HI , 10X , J , ICX , 1HK , 10X, 1HL/ ) 
CD 930 ICTZ=1 , 1 Z 
M - IZQNEt ICTZ, 91 
N = IZONEt ICTZ, 101 
CO 920 J=1 ,N 
CO 910 I =. 1 , M 
II = I * 1 
JJ = J M 

I NO = NODE ( I CTZ , I * J) 

JNO > NODE (ICTZ, II, J) 

KNCI = ' NODEt ICTZ, II, JJ) 
l NO = NODEt ICTZ, I»JJ 1 
WRITEt IOUT, 901) I ELM, I NO, JNO , KNO,LNO 
C 

C PUNCHED OUTPUT IN FORMAT .FOR USE IN PROGRAM SAP 
IF ( IP .EQ. 1) PUNCH 902, I ELM , I NO , JNO ,KNC, LNO 
C 

901 FORMAT (16, 4(6X, I?)} 

902. FORMATt 515) 

. IELM = I ELM + 1 

• 910 CONTINUE 
•920 CONTINUE 
930 CONTINUE 
C . ‘ 

CALL FEMPLTt IZ, INODE, NCOR) 

C 

C * WRITE THE X AND Y COORDINATES OF THE NODE NUMBERS. 

C 

WRITE! I OUT ,94U ) 

940 FORMATt 1H 1 , 4HN0CE , 1 1 X , 1HX , 14 X , 1HY/ ) 

CO 950 1 = 1 ,NCCR. . 

‘ WR ITE ( IOUT ,941 1 INODE, X2t I ) , Y2( I ) 

941 FORMATt 1 H , I 5 ,2 1 5 X,F 1 0.4 1) 

' INODE- = INODE ’+ I 

950 CONTINUE, 

C # 

STOP 
, ENO 


oo 

Cn 


001036 
CCl 04 0 




r 


SUBROUTINE F I NO I LP1 , LP2 , 1 SO} 


U 

c 


SEARCHES PREVIOUSLY CALCULATED DATA TO FIND NODE NUMBERS ASSIGNED TO A 


c 

r 


ZONE BOUNDARY. 

OCOOCfc 

L 


COMMON I CTZ, I ZONE 115, 101 . NODE ( 15 , 25 ,2 5 1 , T EMP ( 25 » 

CC0CC6 

r 


COMMON X2I4002J ,Y214002» 


O 

c 

r 


LP LOCATION OF PT..QN ZONE ICTZ 

0C00C6 

t 


IA = 0 

000007 

r 


IT - 0 

OCOOIO 



CO 5 1=1,25 

000011 



TEMPI I 1 ' = 9999 

000013 

r 

5 

CONTINUE 


L> 

c 

r 


CEFINE CORNOR NOCES ON ZONE' ICTZ 

0C0015 

L 


J1 = IABSI I Z CNE ( ICTZ,LP1)1 

000021 

r 


J2 = IABSI IZ0N6 ( ICTZ , LP 2) ) 

000024 

L 


IIZ = ICTZ 

CC0C25 


10 

IIZ = IIZ - l ■ 

000027 



IF! IIZ J200.200, 110 

CC0030 

r 

110 

CONTINUE 


L 

c 


SEARCH DATA OF ZCNE IIZ 

OOQUJQ 

u 


DO 40 1=1, 7,2 

CC0C32 



11=1+2 , 

0001)34 


• 

IF! 1-7116,15,16 

CC0036 


15 

II = 1 

000037 


16 

CONTINUE 

000037 



K 1 = IABSIIZGNEI III, II) 

000043 

r 


K2 = IABSI IZQNEIIIZ, III) 


l 

c 

r 


COMPARE ICTZ TG IIZ 

0C0047 



IF! J1-K1J 30,20, 3C 

000051 


20 

IFC-J2-K2 140, 21, 40 

000054 


30 

IF(Jl-K2)4G, 31,40 



000056 

CC0Q6C 

000060 

0C0C63 

000063 

000065 

000065 


000066 

000066 

0C0C71 

000072 

000074 

000075 

0C0076 

0C0077 

000101 

0C0IC3 

0C0105 

000110 
0CC11Q 
000112 
OCOl 13 

CC01 16 
000116 
000120 


31 I F { J2-K 1) 40, 21 , 4,Q 
21 CONTINUE 

IA = { 1+11/2 
GO TO 41 ' 

40 CONTINUE 


C 

41 CONTINUE 

If { 14)45,10,45. 

C 

C 

C PUT DESIREC CONTENTS CF IIZ IN TEMPERORY ARRAY 

C 

C 


c 

TEMP ARRAY MUST 

HAVE 

REVERSE .ORDER 

IF 

« • 

c , 

ISD = 1 

CR 2 

AND • IA « 1 

OR 

2 

c 

c 

ISD = 3 

OR 4 

OR - 

AND IA = 3 

OR 

4 


C 

C 

45 CONTINUE 

• MMT = : IZONE< IIZ, 91 ♦ 1 
NNT = I ZONE! IIZ, 10) *■ 1 
MK = MMT 
NK: = NNT 
II = 0 
C 


CO 100 1=1,25 
C 

IF! IA-1146 ,50,46 

46 I F < IA-2 147,60,47 

47 IFt I A- 3)48, 70, 46 
■48 IFI 14-4)100,80,100 

C 

50 CONTINUE 
PK = I 1 • 

NK = 1 

I F ( ISC-2190, 90, 92 
C 

60 CONTINUE 
NK- = I 

IF 1 1 SD-2 ) 91, 91 , 92 


000123 


C 


70 CONTINUE 



000123 



MK = I 

92, SC 

000125 

C 


IF ( I SO-2) 92» 


0C013C 


80 

CONTINUE 


000130 



NK = 1 


0C0131 



NK = I 


000133 

c 


IF! ISC-2192, 

92,91 


c 




000136 


90 

I I = MMT ♦ 1 

- I 

0C0141 



GO TO 93 


000141 


91 

1 1 = NNt + 1 

- I 

000144 


■ 

GO TO 93 


0CC144 


92 

II = I 


000146 


93 

CONTINUE 


000146 



'IF! 11 )200,200,94 

000150 

c 

94 

CONTINUE 

1 

■ 

000150 



TEMP III) = 

NODE! IIZ.Mk.NKI 

000157 


100 

CONTI NUE 


000161 


200 

CONTINUE 


000161 



RETURN 


000162 



END , 



CO 

CO 



000006 

0C0CC6 


000006 

CC0011 

000016 

000022 

000023 

0C0027 


000032 

000034 

000036 

000037 

0C0041 

000042 

CC0044 


GC0G45 

000053 

000061 



C 

c 

c 

c 


c 

c 

c 


. c 


c 


c 

c 


c 


c 


c 

c 

c 


c 

c 

c 


SUBROUTINE F EMPLT 1 1 Z , I NODE jN COR ) 
PLOTS THE FINITE ELEMENT MESH 


COMMON ICTZ, I ZONE 115, 10 1, NODE! 15 , 2 5,2 5) , TEMPI 25) 
COMMON X2I4002) .Y2I4002) 


SCALE DATA 

CALL ASCALEI X2 , 25. *NC0R»1'»20. ) 
CALL ASCALEIY2, 13», NCGR, 1,20 . } 

XSCALE = X2INC0R + 2) 

Y SCALE = Y2I NC OR + 2) 

IFUSCALE .GE. YSCALE ) SF=XSCALE 
IFIYSCALE .GE. XSCALE) SF=YS CALE 


co ioo icri=i,iz 

M = l ZONE ( 1C TZ , 9 ) 

N = I ZONE I ICTZ, 10) 

CO 90 J= 1 , N 
DO 80 1=1, M 

11=1+1 
JJ = J + 1 

CEFIN.E THE 4 NODES OF AN ELEMENT 

. I NO = NODE (ICTZ, I, J) - INODE + 1 
JNO = NODE 1 1 CTZ , 1 1 , J ) - INODE + 1 
KNO = NODE 1 1 CTZ , 1 1 , J J ) - INOCE + 1 
LNO = NODE! ICTZ , I , J J ) t INODE + 1 

I 

DEFINE THE X ANC v' COORDINATES OF THE- 4 NOOES 

: 

XI = X2 { INC) / S F- 
■ XJ = X2(JN0)/SF 
XK = X2IKNC) /SF' 

i 


oo 

VO 


CC0075 

0C0100, 

000102 



00010* 

■C 

XL * 

X2I LNOi /SF 

000106 


VI = 

Y 2 { INOJ/SF 

ocouo 


YJ = 

Y2( JNO/SF 

000112 


VK = 

Y2IKN01/SF 

ccctis 

r 

YL = 

Y2I LNCI/SF 


i. 

c 

/*• 

FLUT 

THE 4 NODES', 

000117 

y 

CALL 

CALPLTIXI ,YIt3) 

000121 


CALL 

CAL PLT ( XJ » YJ > 2 ) 

0C0124 


CALL 

CALPLTI XK,YK,2) 

CC0127 


CALL 

CAL PLT I XLf YLt 2 ) 

000132. 

r 

CALL 

CALPLTI XI i YI ,2 ) 

000135 

u 

80 

CONT I 

NUE 

000142 

90 

CONTINUE 

0CG144 

100 

CONTINUE 

000147 


RETURN 

0CQ150 


END 

• 


VD 

O 
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Table 1. The shape functions. 


x 


8 

E 

i=l 


N. 

i 


x . 
i 


8 

= L 

i=l 


N. y. 
1 x 


n x = - | (i - 5) (l - n) (C + n .+ i) 
,n 2 = j (l -V) (l - n) 
n 3 = | (i +X) (l - n) (5 - n - l) ' 
n 4 - = | (1 + ?) (l- n 2 ) 

n s = | (l + C). (l + n) (5 + n -• i) 

n 6 = | (l - ? 2 ) (l + n) " 
n 7 = j (i - 5 ) (l + n) (-£ + n - l) 

n 8 = | (1 - £) (l - n 2 ) 
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Table 2. Input Data for Example Problem Shown in Figure 5 


EXAMPLE PROBLEM 


3 

lb 

8 ' 

100 

1000 






I 

2 

3 

7 

11 

10- 

9 

6 

3 

. 3 

—3 

4 

5 

8 

13 

12 

-11 

-7 

2 

3 

-11 

-12 

-13 

15 

18 

17 

16 

14 

2 

4 

1 


0. 


0. 






3 


3. 


0. 






5 


5. 


0. 






9 


0. 


3. 






11 


3. 


3. 






13 


3. 


3. 






16 


3. 


7* 






lb 


5m 


7. 
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Table 3. Output for Example Problem Shown in Figure 5 


EXAMPLE PROBLEM 


NO. OF 

ZONES 

3 









ZONE 



ZONE 

NODES 




M 

N 

NO. 

1 

2 

3 

4 

5 

6 

7 

8 


. 

1 

1 

2 

3 

7 

11 

10 

9 

6 

3 

3 

2 

-3 

4 

5 

8 

13 

12 

-11 

-7 

2 

3 

3 

11 - 

12 

-13 

15 

18 

17 

16 

14 

2 

4 


NODE 

X 

Y 

1 

0.000 

0.000 

3 

3.000 

0.000 

5 

5.000 

0.000 

9 

0.000 

3.000 

11 

3.000 

3.000 

13 

5.000 

3.000 

16 

3.000 

7.000 

18 

5.000 

7.000 


(cont ' d . ) " 
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Table 3. Output .for Example Problem Shown in 
F-rgure 5 (continued) . 


ELEMENT 

1 

J 

K 

L 

1000 

100 

101 

105 

104 

1001 

101 

102 

106 

105 

1002 

102 

103 

107 

106 

1003 

104 

105 

109 

108 

1004 

105 

106 

110 

109 

1005 

106 

107 

111 

110 

1006 

108 

109 

113 

112 

1007 

109 

110 

114 

113 

looe 

110 

111 

115 

114 

1009 

103 

116 

118 

107 

101C 

116 

117 

119 

118 

1011 

107 

118 

120 

111 

1012 

118 

119 

121 

120 

1013 

111 

120 

122 

115 

1014 

120 

121 

123 

122 

1015 

115 

122 

125 

124 

1016 

122 

123 

126 

125 

1017 

124 

125 

128 

127 

1018 

125 

126 

129 

128 

1019 

127 

128 

131 

130 

102C 

128 

129 

132 

131 

1021 

130 

131 

134 

133 

1022 

131 

132 

135 

134 


(coht ' d. ) 
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Table 3. Output for Example Shown in 
Figure 5 (concluded) . 


W 00E 

IOC 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 
iL6 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 
127 
•128 

129 

130 

131 

132 

133 

134 

135 


X 

0,0000 

1.0000 

2.0000 

3.QOCO 

0.0000 

1.0000 

2.0000 

3.0000 
0.0000 

1.0000 
2.0000 

3.0000 
0.0000 

1.0000 
2.0000 

3.0000 

4.0000 

5.0000 

4.0000 

5.0000 

4.0000 

5.0000 

4.0000 

5.0000 

3.0000 
4;oogq 

5.0000 
3.0C0Q- 

4.0000 

5.0000 
3.-0000 

4.0000 

5.0000 

3.0000 

4.0000 

5.0000 


y 

0.0000 
0.0000 
0.0000 
0.0000 

1.0000 

1.0000 

1.0000 

1.0000 

2.0000 

2. .00 00 

2.0000 

2.0000 

3.0000 

3.0000 

3.0000 

3.0000 
0.0000 
0.0000 
1 .0000 

1.0000 

2.0000 

2.0000 

3.0000 

3.0000 

4.0000 

4.0000 

4.0000 

5.0000 

5.0000 

5.0000 

6.0000 

6.0000 

6.0000 

7. .00 00 

7.0000 

7.0000 



Table 4. Input for Shear Panel. 


SHcAK HAWEL 


4 

21 • 

17 1 

1000 

I 

2 

3 4 

5 

-11 

-a 

-5 6 

7 

-1 

-10 - 

11 15 

19 

-11 

-12 - 

13 16 

21 

1 

0 . 

0. 


2 

3.0 

0 . 


3 

4.35 

0 . 


4 

4-3727 

.1040 


5 

4.435 

.185 


7 

4. a 

.3 5 


a 

3. 825 

.o25 


y 

3.975 

.975 


li 

2.31 

2.31 


13 

2.475 

2.475 


14 

0 . 

3. 


15 

.825 

3.825 


lo 

.975 

3.975 


17 

0 . 

4. 35 


la 

.1040- 

4.3727 


19. 

.185 

4.435 


21 

-.350 

4.o 



8 

11 

10 

20 

10 

9 

13 

12 

20 

3 

18 

17 

' 14 

10 

20 

20 

-19 

-15 

3 

20 
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Table 5. Input for Bolted Joint Specimen. 

BOLT EU JOINT COMPOSITE STRUCTURE 


15 

66 

41 100 

1000 

23 

29 

30 ' 21 

3 

-30 

22 

7 6 

5 

—7 

-22 - 

30 31 

32 

-9 

-23 - 

32 24 

13 

-30 

41 

57 56 

55 

49 

50 

51 58 

-57 

47 

48 - 

49 -40 

-30 

-51 

-50 - 

49 oO 

61 

-53 

-59 - 

61 43 

-32 

-61 

62 

o 3 44 

34 

-32 

-33 — 

34 25 

15 

-34 

35 

36 2 o 

17 

-36 

37 

38 2 7 

19 

-63 

64 

65 45 

-36 

—•65 

66 

67 

-38 

1 

0 . 

1.5 


3 

0.13 

1..5 


• 4 

0.4290 

1.5 • 


5 

0.5345 

1.5 


6 . 

0.5413 

1. 4635 


7 

0.5625 

1.4325 


8 

0.63 

1.4045 


9 

0.6975 

1.4323 


10 

0.7182 

1. 4635 


11 

0.7255 

1.5 


12 

0.3308 

1.5 


13 

1 . 13 - 

1.5 


15 

1.5 

1.5 


17 

3.3 

1.5 


-19 

4.5 

1.5 . 


22 

0.4292 

1.284 


23 

0.8303 

1.-284 


23 

0 . 

■ 1.0 - 


30 

0.13 - 

1.0 


32 

1-13 

1.0 


34 

1 . 5 - 

1.0 


36 ' 

' 3.3 

1.0 


38 

4.5 

1.0 


41 

0.4292 

0. 7003 


42 

0. 8308 . 

0. 7006 

- 

47 

0 . 

0 . 


49 

0.13 

0 . . 


50 

0.4292 

0.2992 


51 

0.5625 

0.4325 


52 

0.63 

0.4045 


53 

0.6975 

0.4325 


54 

0.7255 

0.5 


55 

0.6975 
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f a) UNGRADED 


(b) GRADED MESH 


Figure 2. Ungraded and Prided Mesh Generated in the Cartesian System. 
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Figure 5. Example of Mesh Generation for a 
Simple Structure. 
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Figure 7 . 


Mesh Generated for Quartersection of Shear Panel. 
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FINITE ELEMENT ANALYSIS OF A COMPOSITE 
BOLTED JOINT SPECIMEN 

By 

Earl A. Thornton 
INTRODUCTION 

With high strength and weight savings/ advanced composite 
materials have become increasingly important in aircraft struc- 
tural design. The full potential for the increase of structural 
efficiency through the use of advanced composites has not yet 
been fully realized because of low efficiencies in mechanical 
joints. The Advanced Composites Design. Guide - (ref. 1). states 
'that weight savings may be reduced by as much--.as 40 percent 
due to such practical constraints. 

. In the use of conventional materials, design methods for 
joints have- evolved ..over a period of time from- data gathered" 
from experimental and. analytical solutions and, in addition,, 
are often based upon rules-of-thumb derived from experience. 

For advanced composites, such data and experience are' relatively 
limited. To partially fill this need, test programs are underway 
at Langley Research Center (LRC) to establish data- on a number 
of mechanical joint designs (ref. 2) 

The purpose of the present study was to provide, analytical 
support for the LRC bolted joint test program. Specific objec- 
tives of the study were to: (1) determine the laminate stress 

distribution in an extra graphite reinforced bolted -joint specimen, 
and (2) compare two methods of modeling bolt transfer loads for 
determination of stress distributions in bolted joints 

This paper will describe the finite element model used to 
represent the bolted joint specimen. The two methods used to 
represent bolt transfer loads will be discussed. Laminate 
membrane force distributions predicted by the finite element 
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analysis will be presented, and force gradients at the bolt holes 
will be discussed. Differences in the results due to the methods 
of representing the bolt loads will. also be ■ discussed. 

BOLTED JOINT SPECIMEN 

The specimen analyzed in this study is the specimen denoted. _ 
as extra graphite reinforced joint specimen number one, reference 
2. The specimen is shown schematically in figure 1 with the 
dimensions used in the analysis. The specimen was fabricated 
from a basic layup of 15 plies reinforced by additional plies 
so that in the thick section where the bolt holes are located 
there are 49 plies. The ply stacking sequences are shown in 
figure 2 with cross-sectional details of the iayup. Reinforcing ■ 
plies increase by 0.1 in. in length per ply over the transition- 
s<=>rvMrm frnm 4Q r>];ieS to 15 plies. 


ANALYTICAL PROCEDURES 
Finite Element Model 

The NASA Stru ctural Analysis (NASTRAN)' computer prograir 
(level 15.5) was used to compute the laminate stress distri- 
butions. in the specimen. The specimen was assumed to be in- 
plane stress and due to symmetry only one— half of the specimen 
was represented with finite elements. The finite element repre- 
sentation is shown in figure 3. The specimen was represented by 
an assemblage of 349 quadrilateral and triangular membrane elements. 
The NASTRAN finite elements used have constant stress throughout 
each element. The mathematical model has -307 grid points and ; 

573 degrees of freedom. Vertical displacements were set to 
zero on the top boundary of the finite element model to repre- 
sent symmetry, and. horizontal displacements at the right edge 
of the finite element model were set to- zero to represent 
clamping in the test fixture. 

In the analytical formulation underlying the present NASTRAN 
elements the element material is assumed homogeneous through the 
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thickness. The element extensional stiffnesses are obtained inter- 
nally in NASTRAN by multiplying the material elasticity matrix by 
the thickness of the element, reference 3. However, the specimen 
in the present study is characterized by several layers of 
material which are assumed homogeneous within the individual 
layers only. Thus for the composite laminate the extensional 
stiffnesses, A^ j , were computed externally using laminated 
plate theory. The stiffnesses were then input to NASTRAN in 
place of the material elasticity matrix-, and the thickness of 
the specimen was everywhere taken as unity. 

The extensional stiffnesses A. . , a 3 x 3 symmetric matrix, 
were computed from reference 4: 


N 

A -^'= £ 
k=l 


-13 




( \ - 


Vi 


in' 


where (Q^j)^ denotes the material elasticity matrix for a 
single layer and (Z^ - denotes -the thickness of- the' 

kth layer. The extensional stiffnesses relate the ih-plane 
membrane forces (N_, N , N . ) to the midplane extensional 

A jL **0 : 

strains (e , e , y ) of- the laminate. Since the extensional 
. y -K-y 

stiffnesses were input to NASTRAN in place of the NASTRAN 

material elasticity matrix, the NASTRAN membrane element 

stresses (o^, -Y xv ) were -the laminate stress resultants 

(N , N , N ) . 
x' y ' xy 

In the present analysis the lamina elastic constants were 
taken as E n =-20 x 10 6 psi, E 22 = 2 x 10 6 psi, G = 0.8 r 
10 6 psi and v^ 2 = 0.3. Each lamina had a thickness of 
0.00542 in. To represent the tapered character of the specimen, 
extensional stiffnesses were computed for the 19 different 
cross-sectional layups. The values of the extensional stiff- 
nesses for the specimen are given in table 1. 
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Bolt Loads 


The specimen was analyzed for loading corresponding to the 
design failure load. This loading, estimated at 21 813 lb was 
assumed to be equally distributed to the three bolts such that 
the total load transmitted to the specimen per bolt was 7271 lb. 

In the finite element model, one-half of this load was applied to 
the center bolt hole and the full value was applied to the lower 
bolt hole. 

Two methods were used to represent the transfer of the bolt 
forces to the finite element model. In the first approach the 
bolt was assumed to have a perfect fit, and the load transfer 
was assumed to take place over 180° of. the bolt hole. .The 
contact force was assigned to vary sinusoidally over this area' 
of contact.. Equilibrium of the bolt was then used to obtain the 
relation: 

2 2Q ' ' 

N = ?T cos 9 (.2) 

where N denotes the contact force per unit arc length, 2Q is 

the total bolt load/ and R -is the radius of the bolt hole. ' The 

angle 8 is measured from a horizontal axis through the hole 

Equation (2) was used to compute equivalent grid point forces 

for each grid point in the contact region (-90° < 0 < 90°). 

/ 

The equivalent grid point forces were computed by integrating 
Equation (2) through an angle of -6° to +6° at each grid ■point 
The equivalent grid point loads are shown in figure 4. 

-In the second approach an imperfect fit was assumed and a 
nonlinear analysis of the bolt transfer loads was made. This - 
analysis, made using the computer program CONTACT -developed in 
reference 5, consists of increasing the bolt load in increments 
and determining the number of grid points in contact and their 
loads at each load increment. The analysis requires as part of 
its input the flexibility matrix for the bolt hole. This flexi- 
bility matrix was obtained from the finite element model by 
applying unit loads at each node of the center bolt hole. The 
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16 x 16 flexibility matrix was computed one column at a time' for 
16 unit load subcases. This matrix was then input to the CONTACT 
program and the bolt transfer forces were computed for several 
load increments . ■ The bolt transfer forces and the region of 
contact for four load increments including the maximum load 
are shown in figure 5 . These forces were computed using an 
initial lack of fit of -0.00287- in. This value, as defined 
in the program, denotes a clearance based upon the radius of 
the hole. 


RESULTS AND DISCUSSION 


The membrane force distributions at the center ' and outside 
bolt holes as predicted by the finite element analysis are 
shown. in figures 6 ' through 8. Shown are plots of the radial ‘ 
force N r , the circumferential force N g , and the in-plane 
shearing force N rQ . versus the angle 0 ■ from the centerline.. 
Predictions based upon the two methods of representing the boll 
transfer loads are compared. 

There is very little, -if any, difference in the membrane 
forces between the center bolt hole and the outside bolt holes, 
Each bolt was assumed to carry the same bolt load and there 
appears to be no interaction effects between holes nor edge 
effects upon the stress distributions in the outside holes. 

The magnitudes and variations of the membrane forces and the 
effects of the two methods of representing the bolt transfer . 
loads can thus be discussed with regard to either hole. 

The largest radial force intensity (fig. 6) occurs, as might 
be expected, on the centerline of the bolt hole. The nonlinear 
bolt loading method predicts the largest radial membrane forces 
with a value of 32 kips/in. compression which is about 23 percent 
higher than the value based upon the cosine bolt loading. The 
largest circumferential membrane force (fig. 7) of 30 kips/in. 
tension occurs at an angle of about 75° from the bolt centerline 
and is also predicted by the nonlinear bolt loading technique. 
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This stress is about 15 percent higher than the value based upon 
the cosine bolt loading. The in-plane membrane shear forces 
(fig. 8) tend to be smaller than the radial or circumferential 
membrane forces. The largest membrane shear force is about 
9 kips/in. and is due to the nonlinear bolt loading. Since 
the in-plane shearing forces tend to be small, the principal 
values (not shown) of the membrane forces correspond in magni- 
tude and location to the maximum radial and circumferential 
membrane forces . 

The distribution of the longitudinal membrane force, N 
along the specimen centerline is shown in figure 9. At x = 0 
the membrane force should be zero since this edge is stress free; 
the small nonzero value is indicative of the error in the finite 
element solution. The membrane force at the left edge of the 
bolt hole (x = 0.5) rises very sharply due to the -indirect 
bearing load of the bolt. On the right side of the bolt hole, 
the force should also be zero since the bolt' is not in contact 
at this point. The finite element solution tends to zero at 
this point. Away from the hole for increasing x, the membrane 
force approaches a uniform value given by the total applied force 
(21 / 816~lb) , divided by the specimen width (3 in.). 

- . Further insight into the results of the finite element 
analyses can be obtained by considering an elasticity solution 
for an isotropic medium. In reference 6, Bickley presents the' 
plane stress- elasticity solution for" a hole in an infinite medium 
loaded by a cosine pressure" distribution over one-half of the 
boundary of the hole.. Closed form solutions for the stress 
components are given in polar coordinates in terms of the 
radius of the hole and Poisson's ratio. Tabulated data of the 
stress components for Poisson's ratio of 0..25 are also presented. 

In figure 10 are shown the membrane force distributions 
predicted by Bickley for an infinite .isotropic medium with a 
hole equal in radius to the bolt hole in the composite specimen 
and loaded by the bolt load used in the finite element analysis. 
The plots are made for r/a =1.2 which corresponds closely to 
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laminated composite material was represented in NASTRAN as a- 
homogeneous material with equivalent extensional stiffness. 
Laminate membrane force distributions were predicted. 

Comparison of the two methods of representing the bolt 
transfer loads showed the two methods were in qualitative 
agreement. The nonlinear analysis estimated membrane forces 
about 20 to 25 percent higher than the linear analysis. Peak 
forces were found to be a radial compressive force on the bolt 
centerline and a circumf erential tensile force of the same magni- 
tude at about 70° from the centerline. In— plane shear forces were 
found to be relatively small. There were little or no interaction 
effects between holes or boundaries of the specimen. Comparison 
of the finite element solution with an isotropic .elasticity 
solution suggests that as a rule these effects will not be 
important for in-plane membrane forces provided offset- distances 
between holes or edges are greater than five hole radii. 
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Table 1 


Section 

- 

An 

A 1 2 

r 

1 

3.01 

0 .•93'2 
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2.94 

0.877 

3 

2.87 

’0. 82'3 

4 

2.65 

0.816 

5 

2.43 

. 0.809 

6 

2 .36 

0.755 

7 

2.28 

■ '0.699 

8 

2.07 

0.693' 
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0.639 

10 

1.92 ' . 
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11 

1.85 

■ 0.529 

12 

1.78 
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13 
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0.419 

14 
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0.364 

15 
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" 0.358 

16 

1.34 ' 

6.303 

17 

1.13 

0.297 

18 

' 1.0'5 

0.242 

19 

1.05 

' 0.242 
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Extensional . stiffnesses 


A. .. x 10“ 6 


(lb/ in. ) 
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- 0.0492 
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0.697 
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0.625 
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0.459 
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0.0492 

0.437 

0.0492 
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Figure 2. • Ply Stacking Sequence for Graphite 
Reinforced Specimen 1. 




MEMBRANE (PLANE STRESS) ELEMENTS 


Figure 3. .Finite Element Representation of Bolted Joint Specimen. 
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Figure 7. Circumferential Membrane Force Distributions at the Bolt Holes 
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(a) INFINITE ISOTROPIC MEDIUM WITH 


COSINE BOLT LOADING 



(b) RADIAL MEMBRANE FORCE DISTRIBUTION FOR r/a ■ 1.2 


Figure 10. 


Elasticity Solution for Infinite 
Me drum with Cosine Bolt Loading. 
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Figure 10 (concluded) . Elasticity Solution for Infinite 

Isotropic Medium with Cosine 
Bolt Loading. 
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